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The  Consolidated  Compartment  Fire  Model  (CCFM)  Computer  Code 
Application  CCFM. VENTS  - Part  III:  Catalog  of  Algorithms  and  Subroutines 


Leonard  Y.  Cooper  and  Glenn  P.  Forney,  Editors 


ABSTRACT 


A project  was  carried  out  at  The  National  Institute  of  Standards  and 
Technology  (NIST)  to  study  the  feasibility  of  developing  a new- generation , 
multi -room,  compartment  fire  model  computer  code,  called  the  Consolidated 
Compartment  Fire  Model  (CCFM)  computer  code.  The  idea  was  that  such  a code 
would  consolidate  past  progress  in  zone -type  compartment  fire  modeling,  and 
allow  readily  for  integration  of  future  advances  with  the  greatest  possible 
flexibility.  Desired  features  of  the  CCFM  would  include:  comprehensive 
documentation,  user-friendliness,  significant  modularity,  numerical 
robustness,  and  versatility  in  the  sense  that  the  code  would  provide  a 
capability  of  analyzing  a particular  compartment  fire  problem  by  using  any  one 
of  a range  of  physical-phenomena-modeling  sophistication,  from  the  most  basic 
to  the  most  comprehensive.  The  project  led  to  the  development  of  a prototype 
multi -room  CCFM  product  called  CCFM. VENTS.  CCFM. VENTS  involves  a model 
formulation  and  code  structure  that  allows  for  the  required  future  CCFM  growth 
flexibility.  It  has  a relatively  sophisticated  and  very  general  room-to-room 
forced  and  unforced  vent  flow  capability.  Finally,  the  CCM. VENTS  code  uses 
the  simplest  possible,  point- source -plume , smoke-filling  fire  physics  in  the 
rooms-of-fire-origin  and  a very  simple  heat  transfer  calculation  there  and  in 
other  spaces. 

This  is  Part  III  of  a four-part  report  which  documents  CCFM. VENTS.  It  is  a 
catalog  of  all  the  modular  algorithms  and  associated  computer  subroutines  used 
to  simulate  the  physical  phenomena  in  CCFM. VENTS.  The  algorithms/subroutines 
are  divided  into  two  categories,  those  which  describe  the  calculations  for 
simulating  the  basic  physical  phenomena  and  those  which  are  used  to  support 
these  calculations.  Tht  latter  catalog  entries  are  called  "utility” 
algorithms/subroutines . 

Each  physical  algorithm  entry  of  the  catalog  includes  a brief  description  of 
the  phenomenon  simulated,  a concise  but  complete  presentation  of  the  calcula- 
tion procedure  used,  identification  of  all  input  and  output  parameters,  and  a 
listing  of  the  subroutine.  The  catalog  entries  have  been  developed  and  are 
presented  as  modular,  stand-alone  products.  The  stand-alone  design  feature 
allows  the  catalog  entries  to  be  used  both  in  CCFM  and  in  any  other  modular, 
zone -type,  compartment  fire  model  computer  code. 

The  other  three  parts  of  this  report  are:  Part  I:  Physical  Basis;  Part  II: 
Software  Reference  Guide;  and  Part  IV:  User  Reference  Guide. 

Keywords:  building  fires;  compartment  fires;  computer  models;  fire  models; 

mathematical  models ; vents ; zone  models . 
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1 . Background 


A project  was  carried  out  at  The  National  Institute  of  Standards  and 
Technology  (NIST)  to  study  the  feasibility  of  developing  a new-generation, 
multi -room,  compartment  fire  model  computer  code,  called  the  Consolidated 
Compartment  Fire  Model  (CCFM)  computer  code.  The  idea  was  that  such  a code 
would  consolidate  past  progress  in  zone -type  compartment  fire  modeling,  and 
allow  readily  for  integration  of  future  advances  with  the  greatest  possible 
flexibility.  Desired  features  of  the  CCFM  would  include:  comprehensive 
documentation,  user-friendliness,  significant  modularity,  numerical 
robustness,  and  versatility  in  the  sense  that  the  code  would  provide  a 
capability  of  analyzing  a particular  compartment  fire  problem  by  using  any  one 
of  a range  of  physical -phenomena -mode ling  sophistication,  from  the  most  basic 
to  the  most  comprehensive.  The  project  led  to  the  development  of  a prototype 
multi -room  CCFM  product  called  CCFM. VENTS.  CCFM. VENTS  involves  a model 
formulation  and  code  structure  that  allows  for  the  required  future  CCFM  growth 
flexibility . 

This  is  the  third  of  a four-part  report  which  documents  all  aspects  of 
CCFM. VENTS. 

In  Part  I [ 1 ] 1 , introductory  remarks  discuss  the  generic  features  of  CCFM,  the 
CCFM  development  process,  and  the  specific  features  of  CCFM. VENTS.  The  main 
objective  of  Part  I is  to  present  a comprehensive  description  and  technical 
basis  of  the  governing  equations  used  in  CCFM. VENTS. 

Part  II  [2]  presents  the  generic  and  CCFM. VENTS -specific  features  of  the 
framework  of  CCFM  software.  This  includes  features  of  both  program  and  data 
structures;  numerical  considerations  used  to  treat  the  solution  to  the  CCFM 
equation  set;  and  a presentation  of  methods  of  how  one  would  use  the  CCFM 
structure  to  generate  stages  of  the  computer  model  beyond  CCFM. VENTS. 

Part  IV  [3]  is  the  CCFM. VENTS  user's  reference  guide. 

The  present  Part  III  is  a catalog  of  all  the  modular  algorithms  and  associated 
modular  computer  subroutines  used  to  simulate  the  physical  phenomena  in 
CCFM. VENTS.  The  purpose  of  this  catalog  is  to  allow  the  interested  user  to 
learn  with  relative  ease  how  the  code  carries  out  any  particular  aspect  of  the 
overall  CCFM. VENTS  simulation.  The  catalog  algorithm/subroutine  entries  are 
also  available  for  use  by  people  who  may  be  interested  in  developing  for  their 
own  particular  needs  a compartment  fire  model  computer  code  other  than  the 
CCFM. 


1 Numbers  in  brackets  always  refer  to  the  references  at  the  end  of  this 
report . 
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2.  Physical  and  Utility  Algorithms/Subroutines  - Modularity  of  the  Catalog 
Entries 

The  CCFM. VENTS  algorithms/subroutines  presented  here  include  two  categories  of 
catalog  entries.  The  first  includes  catalog  entries  which  describe  the 
modeling  of  physical  phenomena.  The  second  category  of  catalog  entry,  those 
which  do  not  involve  the  modeling  of  physical  phenomena  but  which  do  support 
the  calculations  of  the  physical  algorithms,  are  referred  to  as  "utility" 
algorithms.  Documentation  of  utility  algorithms  is  limited  to  a listing  of 
their  appropriately  well -commented  subroutines. 

Each  physical  algorithm/subroutine  entry  of  the  catalog  includes  a brief 
description  of  the  phenomenon  simulated,  a concise  presentation  of  the 
calculation  procedure  used,  identification  of  all  input  and  output  parameters, 
and  an  actual  listing  of  the  computer  subroutine.  The  catalog  provides  the 
bridge  between  the  governing  equations  of  CCFM. VENTS  and  their  technical 
basis,  as  presented  in  detail  in  Part  I,  and  the  implementation  of  these 
mathematical  submodels  into  the  software  of  CCFM. 

The  entries  of  the  catalog  have  been  developed  and  are  presented  as  stand- 
alone products.  The  stand-alone  design  feature  allows  the  catalog  entries  to 
be  used  readily  both  in  CCFM  and  in  any  other  modular,  zone -type,  compartment 
fire  model  computer  code. 

Implementing  the  CCFM  concept  would  lead  to  a process  of  development  which 
would  consolidate  past  progress  in  zone -type  compartment  fire  modeling,  and 
which  would  allow  readily  for  integration  of  future  advances  with  the  greatest 
possible  flexibility.  It  is  envisaged  that  the  modularity  and  stand-alone 
features  of  the  catalog  entries  would  help  to  provide  this  flexibility. 

Indeed,  these  features  would  be  critical  in  insuring  the  required  orderly  and 
efficient  future  growth  of  and  improvement  to  existing  stages  of  the  CCFM. 

For  example,  while  some  of  the  algorithms/subroutines  presented  here  and  used 
in  CCFM. VENTS  involve  very  simple  mathematical  models  of  the  physical 
phenomena  simulated,  there  is  clear  opportunity  for  improvement  by  the  use  of 
more  sophisticated  modeling  ideas.  The  modular  design  of  both  the  original 
and  the  new  catalog  entry  and  subroutine  software  would  then  allow  easily  for 
the  replacement  of  the  old  product  with  the  new,  or  for  the  option  of  using 
the  product  of  either  greater  or  less  sophistication,  say,  for  comparing  the 
supposed  benefit  of  using  one  or  the  other.  Examples  of  this  are  the 
presently  used  plume  algorithm,  PLUGO,  which  implements  a very  simple  point- 
source  fire  plume  model,  and  the  presently  used  vent-flow  distribution  model, 
FL0G02 , which  extends  and  implements  a very  simple,  well- tested  rule  for 
distributing  vent  flows  and  calculating  heat  transfer  to  the  upper  and  lower 
layer  and  to  the  surfaces  of  the  receiving  room.  More  sophisticated  alterna- 
tives to  both  of  these  algorithms  are  now  under  development. 

The  use  of  well -documented  modular  algorithms/subroutines  in  CCFM. VENTS  would 
also  be  critical  for  successful  future  CCFM  team  development  efforts.  Indeed, 
by  constructing  CCFM  products  with  modular  components  having  well-defined 
inputs  and  outputs,  a development- team  effort  that  involves  significant  inter- 
institution participation  would  be  possible. 
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3.  The  Format  of  the  Catalog  Entries 


The  development,  technology  transfer,  and  use  of  CCFM  is  enhanced  by  es- 
tablishing and  implementing  guidelines  for  a uniform  format  of  algorithm/sub- 
routine documentation.  In  this  regard  a prototype  format  has  been  developed 
and  used  in  all  CCFM. VENTS  algorithm/subroutine  catalog  entries  presented 
here.  It  is  hoped  that  this  first  CCFM  catalog  will  provide  a useful  test  of 
the  proposed  format,  and  that  this  will  lead  to  effective,  enriched,  future 
editions  of  the  CCFM  catalog. 

The  adopted  algorithm  format  includes  the  following  elements: 


TITLE 


Should  indicate  the  main  purpose  of  the 
algorithm/subroutine . 


DESCRIPTION  - General  description  of  the  algorithm. 

OUTPUT  - List  of  output  variables , including  definitions 

and  units . 


INPUT 


List  of  input  variables,  including  definitions 
and  units . 


CALCULATIONS 


SUBROUTINES 

USED 


Concise  description  of  the  rules  for  obtaining 
the  output  variables  from  the  input  variables . 
This  would  include  or  refer  explicitly  to  all 
equations  required  in  the  calculation.  If  other 
algorithms/subroutines  are  required,  then  these 
should  be  readily  available  and  referenced. 

Listing  of  or  explicit  reference  to  each 
algorithm/subroutine  used  to  carry  out  the 
calculations . 


REFERENCES 

SUBROUTINE 

VARIABLES 

PREPARED  BY 

SUBROUTINE 


A list  of  references. 

A cross-reference  of  all  nomenclature  (including 
units)  introduced  in  the  above  sections  to  the 
nomenclature  used  in  the  FORTRAN  computer 
subroutine . 

Names  of  those  who  prepared  the  algorithm/sub- 
routine and  date  of  preparation. 

Listing  of  the  computer  subroutine.  This  would 
be  well -commented  and  would  include  a summary  of 
the  purpose  of  the  subroutine  and  definitions 
(including  units)  of  its  input  and  output 
variables . 
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All  algorithm/subroutines  used  in  CCFM. VENTS  are  presented  here  in  Appendices 
A and  B.  Appendix  A is  a compilation  of  the  physical  algorithm  catalog 
entries.  All  of  these  are  presented  in  the  above  format.  Appendix  B is  a 
compilation  of  the  supporting  utility  catalog  entries. 


4.  Tabulation  of  the  Catalog  Entries 

The  entries  of  each  appendix  are  compiled  in  alphabetical  order  according  to 
algorithm/subroutine  name,  with  independent  pagination  for  each  entry.  This 
presentation  highlights  the  modular  and  stand-alone  characteristic  of  the 
entries  and  the  flexibility  of  the  overall  CCFM.  Thus,  at  any  time,  new  or 
improved  entries  can  be  sorted  into  or  used  to  replace  existing  entries, 
thereby  leading  to  an  expanded  and  improved  CCFM  catalog.  Also,  as  ap- 
propriate, the  catalog  entries  can  be  lifted  and  used  elsewhere  to  construct 
and/or  enhance  two -layer  zone -type  compartment  fire  models  other  than  CCFM. 

Listings  of  the  catalog  entries  in  both  Appendices  are  presented  in  the  Table 
of  Contents . 
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APPENDIX  A - Physical  Algor i thus /Sub routines 


C0MVL1 


Setup  for  Calculation  of  the  Flow  Through  Natural  Vents  in  the 
Wall  Segment  Common  to  Two  Rooms 


DESCRIPTION 

Consider  an  instant  of  time  during  the  simulation  of  a compartment  fire 
environment.  This  algorithm  will  be  useful  in  setting  up  a procedure  for 
calculating  the  instantaneous  flow  rate  through  one  or  more  natural  vents  of 
arbitrary  size  and  elevation  located  in  the  wall  segment  common  to  any  two 
particular  rooms  of  the  compartment  [1]. 

Designate  the  two  rooms  as  rooms  1 and  2 and  refer  to  the  Figure  1.  In 
general,  within  the  elevation  interval  of  the  common  wall  segment  the  room  1- 
to-room  2 pressure  difference  can  be  described  by  at  most  three  piecewise- 
continuous,  linear  functions  of  elevation.  These  will  have  endpoints  at  the 
bottom  and  top  of  the  common  wall  segment  and  at  the  elevations  of  the 
interfaces  which  separate  the  upper  and  lower  layers  of  the  two  rooms  (provi- 
ded these  interfaces  also  lie  in  the  elevation  interval  of  the  wall  segment) . 
It  is  possible  for  each  of  the  pieces  of  the  piecewise  linear  functions  to 
change  sign  somewhere  within  its  operable  elevation  range.  The  flow  calcula- 
tion for  an  arbitrary  vent  in  the  wall  segment  can  be  carried  out  conveniently 
once  the  elevations  (no  more  than  four)  and  corresponding  room  1 and  2 
absolute  pressures  and  pressure  differences  at  the  endpoints  of  these  pieces 
have  been  identified  explicitly.  (The  abso1  te  pressures  will  be  important  in 
cross-vent  flow  calculations  only  in  relatively  unusual  circumstances  where 
compressiblility  effects  are  significant.)  Also  required  for  the  calculation 
are  the  elevations  (no  more  than  three  in  addition  to  the  above)  where  the 
pieces  change  sign  (i.e.,  elevations  where  the  room  1-to-room  2 pressure 
difference  is  zero) . 

Pressure  distributions  in  each  room  and  possible  cross -wall  pressure - 
difference  distributions  are  depicted  in  Figure  2. 

This  algorithm  calculates  the  total  number  cf  these  elevations  (no  more  than 
seven  unique  elevations) , NELEV ; the  values  :he  elevations  relative  to  a 
datum  elevation  (presented  in  an  array,  incr«_  ^.ing  uniformly  with  elevation), 
yN ; and  corresponding  rrays  of  absolute  pressure  in  rooms  1 and  2,  px  N and 
p2  N , respectively,  the  room  1-to-room  2 pressure  differences,  Apx  2 N. 


OUTPUT 


E L E V 


Number  of  unique  elevations,  < 7,  which  identify  1) 
endpoints  of  the  continuous,  linear  pieces  of  the  room 
1-to-room  2 pressure  difference  vs  elevation  and/or  2) 
elevations  where  the  room  1-to-room  2 pressure 
difference  is  zero,  where  all  these  elevations  are  at 
or  between  the  vertical  extremes  of  the  wall  segment 
common  to  rooms  1 and  2 . 


C0MWL1 


1 


Pi.h  [p2>N] ; n - i to  nwelev 


Absolute  hydrostatic  pressure  in  room  1 [room  2]  at 
elevation  yN . [Pa  - kg/(m-s2)] 


yN;  N « 1 to  Nwelev ; and  yN  < yN+1 

Values  of  the  elevations  identified  under  the  output 
variable  NWELEV  above  the  datum  elevation,  [m] 


aPi,2,n  : N = 1 to  nwelev 

Pressure  in  room  1 - pressure  in  room  2 at  elevation 
yN  . If  Apx  2 N >0,  *=0,  or<0  then  flow  through  a 
vent  at  this  elevation  will  be  from  room  1 to  room  2 
zero,  or  from  room  2 to  room  1,  respectively. 

[Pa  = kg/(m-s2)] 


INPUT 

Pdatum 


Datum  absolute  pressure.  [Pa  ■=  kg/(m-s2)] 


Pfloor  , i > i ■=  1 or  2 

Hydrostatic  pressure  at  floor  in  room  1 [room  2]  above 
the  datum  pressure.  [Pa  — kg/(m-s2)] 


yFLooR.i  tycEiL.il  Flayer, i);  i ■=  1 or  2 

Elevation  of  the  floor  [ceiling]  {layer  interface}  in 
room  i above  the  datum  elevation,  [m] 

Pt , i [Pu. i 1 ! i “ 1 °r  2 

Density  of  lower  [upper]  layer  in  room  i.  [kg/m3] 


CALCULATIONS 

1.  Set  yMIN  to  maximum  of  two  floor  elevations  and  yMAX  to  minimum  of  two 
ceiling  elevations. 


C0MWL1  - 2 


2. 


Store  into  the  elevation  list,  yN , the  floor,  ceiling  and  layer 
elevations  of  the  two  rooms  that  are  between  yMIN  and  yMAX . Use  the 
algorithm/subroutine  DELP  to  calculate  the  absolute  pressure  in  each 
room  and  the  difference  in  room  pressures  at  each  of  these  elevations 
and  store  the  results  in  the  arrays  px  N , p2  N , and  Apx  2 N . 

3.  Sort  the  elevation/pressure/pressure-difference  list  according  to 
elevation  and  remove  any  duplicate  entries. 

4.  Identify  and  add  any  additional  neutral  planes  (elevations  where 
pressure  difference  is  zero)  to  the  elevation/pressure/pressure- 
difference  list. 

5.  Re-sort  list,  remove  duplicate  entries,  and  set  NWELEV  to  number  of 
elevation  entries  in  list.  Return  NWELEV , px _ N , P2 , N . yN , APi , 2 , n to 
calling  routine. 


SUBROUTINES  USED 

DELP  [2] 

HSORT  [3,4] 
RMVDUP  [5] 
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SUBROUTINE  VARIABLES 


All  nomenclature  in  the  subroutine  is  identical  to  the  nomenclature  used  above 
except  for: 


^HELEV 

' 

NELEV 

Pfloor , 1 . Pfloor , 2 

- 

PL0R(1) , 

PL0R(2)  [Pa  = kg/(m- s2 ) 

Pi  , N [P2  , N ] 

- 

PAB1(N) , 

PAB2 (N)  [Pa  = kg/(m-s2) 

yN 

- 

YELEV(N) 

[m] 

YcEIL , 1 > yCEIL  , 2 

- 

YCEIL(l) , 

, YCEIL(2)  [m] 

yFLOOR , 1 > yp  L 0 OR , 2 

- 

YFLOR(l) , 

, YFLOR ( 2 ) [m] 

yLAYER , 1 ’ yLAYER , 2 

- 

YLAY(l) , 

YLAY(2)  [m] 

APl , 2 , N 

- 

DP1M2 (N) 

[Pa  = kg/(m-s2)] 

PL  , 1 ’ Ph  , 2 

- 

DENL(l) , 

DENL(2)  [kg/m3] 

Pu  , 1 > PU  , 2 

- 

DENU(l) , 

DENU(2)  [kg/m3] 

PREPARED  BY 

Leonard  Y.  Cooper  and  Glenn  P.  Forney 
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Figure-  1.  The  fire -generated  environment  in  two  adjacent  rooms 
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Figure  2. 


Possible  distributions  of  Apx  2 for  a common-wall  segment. 
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SUBROUTINE 


SUBROUTINE  COMWLl( 

I YFLOR , YCEIL , YLAY , DENU , DENL , PFLOR , 

I PDATM, 

0 YELEV,  DP1M2 , PAB1 , PAB2 , NELEV 
+ ) 

IMPLICIT  DOUBLE  PRECISION  (A-H.O-Z) 

C*BEG 

C 


C*** 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


COMWL1  GENERIC  - SETUP  FOR  CALCULATION  OF  THE  FLOW  THROUGH  VENTS 
IN  THE  WALL  SEGMENT  COMMON  TO  TWO  ROOMS.  THIS  ROUTINE 
CALCULATES  ROOM  PRESSURES  AND  ROOM-TO-ROOM  PRESSURE 
DIFFERENCES  AT  CERTAIN  ELEVATIONS  ALONG  COMMON  WALLS  OF 
ADJACENT  ROOMS.  THESE  ELEVATIONS  ARE: 

1.  YMIN  - MAXIMUM  OF  TWO  ROOMS  FLOOR  ELEVATIONS 

2.  YMAX  - MINIMUM  OF  TWO  ROOMS  CEILING  ELEVATIONS 

3.  LAYER  ELEVATIONS  BETWEEN  YMIN  AND  YMAX 

4.  NEUTRAL  PLANES  ELEVATIONS  BETWEEN  YMIN  AND  YMAX 

THERE  CAN  BE  UP  TO  7 ELEVATIONS  OF  INTEREST  (1  FLOOR,  1 CEILING, 

2 LAYER  AND  3 NEUTRAL  PLANE  ELEVATIONS) . 


C***  SUBROUTINE  VARIABLES 
C INPUT 
C 

C YFLOR  - ELEVATION  OF  FLOOR  IN  EACH  ROOM  [M] 

C YCEIL  - ELEVATION  OF  CEILING  IN  EACH  ROOM  [M] 

C YI AY  ELEVATION  OF  LAYER  IN  EACH  ROOM  [M] 

C DJNU  - UPPER  LAYER  DENSITY  IN  EACH  ROOM  [KG/M**3] 

C DENL  - LOWER  LAYER  DENSITY  IN  EACH  ROOM  [KG/M**3] 

C PFLOR  - PRESSURE  AT  FLOOR  OF  EACH  ROOM  [KG/(M*S**2)  - PA] 

C PDATM  - ABSOLUTE  PRESSURE  AT  REFERENCE  ELEVATION  [PA] 

C 

C OUTPUT 

C 

C NELEV  - NUMBER  OF  ELEVATIONS  (DESCRIBED  IN  1 -->  4 ABOVE) 

C YELEV  - ELEVATIONS  OF  ELEVATIONS  (DESCRIBED  IN  1 -->4  ABOVE)  [M] 

C DP1M2  - PRESSURE  DIFFERENCES  AT  ELEVATIONS  IN  YELEV  [PA] 

C PAB1  - ARRAY  OF  ABSOLUTE  PRESSURES  AT  ELEVATIONS  OF  INTEREST  [PA] 

C PAB2  - ARRAY  OF  ABSOLUTE  PRESSURES  AT  ELEVATIONS  OF  INTEREST  [PA] 

C 

C***  ROUTINES  CALLED : 

C DELP , HSORT , RMVDUP 

C 

C*END 


DIMENSION  PAB1 (*) , PAB2(*) 

DIMENSION  YFLOR (*) , YCEIL(*) , YLAY(*) , DENU(*) , DENL(*) 
DIMENSION  PFLOR (*) , YELEV (*) , DP1M2(*) 
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YMIN  «=  MAX(YFL0R(1)  ,YFL0R(2)) 

YMAX  = MIN(YCEIL(1) ,YCEIL(2)) 

YELEV(l)  - YFLOR(l) 

YELEV ( 2 ) - YFLOR(2) 

YELEV ( 3 ) - YLAY(l) 

YELEV (4)  - YLAY(2) 

YELEV ( 5 ) - YCEIL(l) 

YELEV ( 6 ) - YCEIL(2) 

C 

NELEV  = 0 
DO  10  I = 1,  6 
C 

c***  COPY  ONLY  THOSE  VALUES  OF  Y(I)  BETWEEN  YMIN  AND  YMAX 
C 

IF (YMIN . LE . YELEV ( I ) . AND .YELEV(I) . LE .YMAX )THEN 
NELEV  = NELEV  + 1 
YELEV (NELEV)  = YELEV ( I ) 

CALL  DELP( 

I YELEV (NELEV) , YFLOR , YLAY , DENL , DENU , PFLOR , PDATM , 

0 DP1M2 (NELEV) , PABI (NELEV) , PAB2 (NELEV) ) 

ENDIF 
10  CONTINUE 
C 

C***  SORT  THE  LIST  OF  HEIGHTS  AND  REMOVE  ANY  DUPLICATES  THAT  ARE  FOUND 
C 

CALL  HSORT (YELEV , DP1M2 , PABI , PAB2 , NELEV) 

CALL  RMVDUP (YELEV , DP1M2 , PABI , PAB2 , NELEV) 

N - NELEV 
DO  20  I - 2,  NELEV 
C 

C***  LOOK  FOR  NEUTRAL  PLANES. 

C 

IF(DP1M2(I)*DP1M2(I-1) .LT.O. )THEN 
C 

C***  DP1M2  HAS  A CHANGE  OF  SIGN,  SO  INTERPOLATE  TO  FIND  THE 
C HEIGHT  WHERE  DP1M2  IS  ZERO. 

C 

C 

C SOLVE  FOR  YINTRP 

C 

C YINTRP  - YELEV ( I - 1 ) 0.  - DP1M2(I-1) 

C - 

C YELEV ( I ) - YELEV ( I - 1 ) DP1M2(I)  - DP1M2(I-1) 

C 

Y2  - YELEV ( I ) 

Y1  - YELEV ( I - 1 ) 

P2  - DP1M2 (I) 

PI  - DP1M2(I-1) 

YINTRP  - (Y2*P1-Y1*P2)/(P1-P2) 

C 

C***  CALL  DELP  TO  CALCULATE  ABSOLUTE  PRESSURES  AT  NEUTRAL  PLANE 
C 
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CALL  DELP( 

I YINTRP , YFLOR , YLAY , DENL , DENU , PFLOR , PDATM , 

0 DUMMY , PAB1 (N+l ) , PAB2 (N+l ) ) 

YELEV(N+1)  - YINTRP 
DP1M2 (N+l)  - 0.0D0 
N « N + 1 
END  IF 
20  CONTINUE 
C 

C***  RE-SORT  LIST  AND  REMOVE  DUPLICATES 
C 

I F ( NELEV . NE . N ) THEN 
NELEV  - N 

CALL  HSORT(YELEV , DP1M2 , PAB1 , PAB2 , NELEV) 

CALL  RMVDUP ( YELEV , DP1M2 , PAB1 , PAB2 , NELEV) 
ENDIF 
RETURN 
END 
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DELP 


Calculation  of  the  Absolute  Hydrostatic  Pressures  at  a Specified 
Elevation  in  Each  of  Two  Adjacent  Rooms  and  the  Presssure 
Difference 


DESCRIPTION 

Consider  an  instant  of  time  during  the  simulation  of  a compartment  fire 
environment.  This  algorithm  calculates  the  absolute  hydrostatic  pressures  at 
a specified  elevation  in  each  of  two  adjacent  rooms  of  the  compartment  and  the 
difference  of  these  pressures  [1] . These  results  are  required  in  order  to 
characterize  the  instantaneous  rate  of  flow  across  a vent  in  a wall  segment 
common  to  the  two  rooms . 

Designate  the  two  rooms  as  rooms  1 and  2.  Using  the  pressure  algorithm  of  [2] 
this  algorithm  calculates  the  absolute  hydrostatic  pressures  in  the  rooms,  px 
and  p2 , respectively,  at  a specified  elevation  and  then  their  difference, 

Apx  2 = Pi  - p2 • Depending  on  whether  Apx  2 is  positive  or  negative,  the 
flow  through  a vent  near  the  elevation  in  question  will  be  from  room  1 to  room 
2,  or  from  room  2 to  room  1,  respectively. 

OUTPUT 


Pi  [P2] 


Absolute  hydrostatic  pressure  in  room  1 [room  2]  at 
elevation  y [Pa  = kg/(m-s2)]. 


APl,2 


Pi  - P2  tPa  “ kg/(m-s2) ] . 


INPUT 


Pdatum 


Datum  absolute  pressure  [Pa  *=  kg/(m-s2)]. 

Pfloor , i > i = 1 or  2 

Pressure  at  floor  in  room  1 [room  2]  above  datum 
pressure  [Pa  = kg/(m-s2)]. 


DELP  - 1 


Elevation  above  the  datum  elevation  where  Apx  2 is  to 
be  calculated  [m] . 


^floor  , i Flayer, il»  i = 1 or  2 

Elevation  of  the  floor  [layer  interface]  in  room  i 
above  the  datum  elevation  [m] . 


Pl, i [Pn.il:  i = 1 or  2 

Density  of  lower  [upper]  layer  in  room  i [kg/m3]. 


CALCULATIONS 


For  rooms  i = 1 , 2 calculate  the  pressure  and  then  the  pressure  difference  at 
height  y above  the  floor. 


6Pi (y) 


Pt.  iS<y  - yFLOOR.i);  /floors  ^ y * yLAYER.i 

Pl , i S(YlaYER , i ' yFLOOR.i)  ' Pu , i S(y  ' yL AYER , i ) » 

yLAYER , i < y 


Pi Cy)  PFLOOR.i  + 5Pi (y)  + PdATUM 

APi , 2 s Pi(y)  - p2(y)  = (pfloor.i  - pf l o o r , 2 ) + [fipi(y)  - 5P2(y)l 


where 


(1) 

(2) 

(3) 


g = 9.8  m/s2 

is  the  acceleration  of  gravity. 


REFERENCES 

[1]  Cooper,  L.Y.  and  Forney,  G.P.,  Section  3.3  of  Consolidated  Compartment 
Fire  Model  (CCFM)  Computer  Code  Application  CCFM. VENTS  - Part  I: 
Physical  Basis,  NISTIR  90-4342,  National  Institute  of  Standards  and 
Technology,  Gaithersburg  MD,  1990. 

[2]  Tanaka,  T.,  A Model  of  Multiroom  Fire  Spread,  NBSIR  83-2718,  U.S. 
National  Bureau  of  Standards,  Gaithersburg,  MD,  1978. 
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SUBROUTINE  VARIABLES 


Pdatum 

PDATM  [Pa  - kg/(m-s2)] 

Pi.  P2 

PAB1 , PAB2  [Pa  - kg/(m-s2)] 

y 

Y [m] 

yFLOOR , 1 1 yFLOOR , 2 

YFLOR(l) , YFLOR ( 2 ) [m] 

yL AYER . 1 > yLAYER  , 2 

YLAY(l),  YLAY(2)  [m] 

APl,2 

DP  [Pa  = kg/(m-s2)] 

5px , Sp2 

PROOM(l) , PROOM ( 2 ) [Pa  = kg/(m-s2) 

Pl  , 1 . Pl.Z 

DENL(l) , DENL(2)  [kg/m3] 

P\i  , 1 > P\i  , 2 

DENU(l) , DENU(2)  [kg/m3] 

PREPARED  BY 

Leonard  Y.  Cooper  and  Glenn  P. 
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SUBROUTINE 


SUBROUTINE  DELP( 

I Y , YFLOR , YLAY , DENL , DENU , PFLOR , PDATM , 
0 DP,  PAB1 , PAB2) 

IMPLICIT  DOUBLE  PRECISION  (A-H.O-Z) 

C*BEG 


C*** 

C 

C 

C 

C 

C 

c*** 

c 


DELP  GENERIC  - CALCULATION  OF  THE  ABSOLUTE  HYDROSTATIC 

PRESSURES  AT  A SPECIFIED  ELEVATION  IN  EACH  OF  TWO  ADJACENT 
ROOMS  AND  THE  PRESSURE  DIFFERENCE.  THE  BASIC  CALCULATION 
INVOLVES  A DETERMINATION  AND  DIFFERENCING  OF  HYDROSTATIC 
PRESSURES  ABOVE  A SPECIFIED  DATUM  PRESSURE. 

SUBROUTINE  ARGUMENTS 


C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 

C*END 


INPUT 


Y - HEIGHT  ABOVE  DATUM  ELEVATION  WHERE  PRESSURE  DIFFERENCE  IS 
TO  BE  CALCULATED  [M] 

YFLOR  - HEIGHT  OF  FLOOR  IN  EACH  ROOM  ABOVE  DATUM  ELEVATION 
[M] 

YLAY  - HEIGHT  OF  LAYER  IN  EACH  ROOM  ABOVE  DATUM  ELEVATION  [M] 
DENL  - LOWER  LAYER  DENSITY  IN  EACH  ROOM  [KG/M**3] 

DENU  - UPPER  LAYER  DENSITY  IN  EACH  ROOM  [KG/M**3] 

PFLOR  - PRESSURE  AT  BASE  OF  EACH  ROOM  ABOVE  DATUM  PRESSURE 
[KG/(M*S**2)  -=  PASCAL] 

PDATM  - ABSOLUTE  DATUM  PRESSURE  [KG/(M*S**2)  *=  PASCAL] 


OUTPUT 


DP 

PAB1 


PAB2 


CHANGE  IN  PRESSURE  BETWEEN  TWO  ROOMS  [KG/(M*S**2)  « PASCAL] 
ABSOLUTE  PRESSURE  AT  ELEVATION  Y IN  ROOM  1 [KG/(M*S**2)  - 
PASCAL] 

ABSOLUTE  PRESSURE  AT  ELEVATION  Y IN  ROOM  2 [KG/(M*S**2)  *= 
PASCAL] 


DIMENSION  PROOM(2) 

DIMENSION  DENL(*) , DENU(*) , YFLOR(*) , YLAY(*) , PFLOR(*) 
PARAMETER  (G-9.80D0) 

DP  - 0.0D0 
DO  10  I - 1,  2 

PROOM(I)  - 0.0D0 

IF(YFLOR(I) . LE . Y . AND . Y . LE . YLAY ( I ) ) THEN 
C 

C***  THE  HEIGHT,  Y,  IS  IN  THE  LOWER  LAYER 
C 

PROOM(I)  - - (Y  - YFLOR(I) )*DENL(I)*G 

C 

ELSE  I F ( Y . GT . YLAY ( I ) ) THEN 
C 

C***  THE  HEIGHT,  Y,  IS  IN  THE  UPPER  LAYER 
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PROOM(I)  - - (YLAY(I) -YFLOR(I) )*DENL(I)*G 

+ - (Y-YLAY(I) )*DENU(I)*G 

END  IF 
10  CONTINUE 
C 

C***  CHANGE  IN  PRESSURE  IS  DIFFERENCE  IN  PRESSURES  IN  TWO  ROOMS 
C 

DPI  - PFLOR(l)  + PROOM(l) 

DP2  = PFLOR(2)  + PROOM ( 2 ) 

PAB1  -=  PDATM  + DPI 
PAB2  = PDATM  + DP2 

C IF(ABS (DPI -DP2) . LT .MAX(1 . 0D- 3 , 1 . OD- 3*(ABS (DP1)+ABS (DP2) ) ) )THEN 
C DP  = O.ODO 

C ELSE 

DP  - DP1-DP2 
C ENDIF 

RETURN 
END 
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FANRES 


Calculate  the  Flow  Through  a Fan/Duct  Forced  Ventilation  System 


DESCRIPTION 

FANRES  is  a model  for  predicting  flow  through  a simple  fan/duct  forced 
ventilation  system  connecting  any  two  rooms  of  a multi-room  facility,  or 
connecting  a room  and  the  outside  environment  [1].  It  is  based  on  an 
extension  to  the  ideas  presented  in  [2].  The  FANRES  system  consists  of  a fan 
and  duct  with  a single  inlet  and  outlet.  The  model  predicts  flow  through  the 
system  under  arbitrary  end-point  conditions.  For  example,  if  an  appropriate 
fan  curve  is  available,  FANRES  can  provide  an  estimate  of  flow  rate  even  when 
the  room  environments  at  the  end-points  lead  to  flow  through  the  system  in  the 
direction  opposite  to  normal  fan  operation. 

Heat  transfer  is  neglected  and  the  fan/duct  system  is  assumed  to  be  adiabatic. 

Refer  to  Figure  1.  Designate  the  two  rooms  of  interest,  or  the  room  and  the 
outside  environment  as  rooms  1 and  2.  Define  room  1 [room  2]  as  the  room 
connected  to  the  system  inlet  [outlet].  Here,  the  system  inlet  [outlet]  is 
the  end-point  of  the  fan/duct  flow  system  which  corresponds  to  the  flow  inlet 
[outlet]  under  normal  fan  operating  conditions. 

Characterize  the  elevations  of  the  inlet  and  outlet  ducts  relative  to  a datum 
elevation  by  the  elevation  of  their  midpoints,  yE y s t e m i anc*  ysYSTEM  2> 
respectively.  Designate  also  the  hydrostatic  pressures  in  the  rooms  at  these 
elevations,  relative  to  a datum  pressure,  as  6pSYSTEM  i anc*  ^Psystem  2> 
respectively.  Then,  it  is  assumed  that  the  system  elevation  and  pressure 
differences,  AySYSTEM  and  ApSYSTEM . respectively,  are  specified,  where 


SYSTEM 

^SYSTEM , 1 

^SYSTEM , 2 

(1) 

^P  SYSTEM 

= ^PSYSTEM , 1 

<5PSYSTEM,2 

(2) 

The  FANRES  algorithm  requires  specified  characteristics  of  the  fan/duct 
system.  These  include  the  fan  operating  curve  and  the  effective  flow 
resistance  of  the  ducting. 

The  Fan  Operating  Curve 

The  fan  curve  provides  the  volume  flow  rate  through  the  fan,  VF  as  a function 
of  cross -fan  pressure  difference,  ApF , i.e., 

VF  = f(ApF)  (3) 

where 

ApF  = <5pF  , 2 " ^Pf  , l 
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and  5pF  x and  5pF  2 are  the  pressures  that  would  be  measured  at  the  fan  inlet 
and  outlet  sections,  respectively.  A sketch  of  a generic  fan  operating  curve 
is  presented  in  Figure  2.  As  indicated  in  the  figure,  a fan  curve  used  in 
FANRES  must  be  such  that  ApF  is  monotonically  decreasing  with  increasing  VF . 

In  terms  of  the  fan  curve,  forward  [backward]  flow,  from  room  1 to  room  2 
[room  2 to  room  1],  is  associated  with  positive  [negative]  values  of  the 
volume  flow-rate,  VF  , and  corresponding  positive  [negative]  values  of  the  mass 
flow- rate,  MF  . 

For  the  purpose  of  the  present  calculations,  the  fan  curve  of  Eq.  (3)  is 
approximated  by  linear  interpolation  between  and  linear  extrapolation  beyond 
two  or  more  pairs  of  specified  values,  [VF  j,  ApF  j ] , J ■=  1 to  JTERM  > 1, 
which  identify  [VF , ApF ] points  on  the  curve.  If  only  one  point  is  specified, 
i.e.,  JTERM  = 1 with  specified  Vp  1 and  ApF  x,  then  FANRES  models  the  fan 
curve  as  a line  of  constant  VF  = VF  ^ . In  the  latter  case,  the  fan/duct 
system  simulated  is  assumed  to  be  one  which  always  delivers  the  specified 
volume  flow-rate. 


Flow  Resistance  of  the  Duct 


Since  the  environments  of  the  two  rooms  1 and  2 are  generally  different,  the 
density  of  the  flow  through  the  system  depends  on  flow  direction,  i.e.,  which 
of  the  two  rooms  is  the  source  of  the  flow.  Flow  resistance  through  the  duct, 
to  be  defined  below,  would  also  be  dependent  generally  on  flow  direction. 
Define  the  density  of  the  flow  and  the  resistance  of  the  duct  as  pDUCT  1 and 
Rx  , respectively,  if  the  flow  is  forward  flow,  from  room  1,  and  as  pDUCT  2 anc* 
R2 , respectively,  if  the  flow  is  backward  flow,  from  room  2.  The  pDUCT  i and 
Ri  are  assumed  to  be  specified. 

From  the  above-specified  parameters  define  and  evaluate  the  adjusted  pressure 
differences,  ApADJ  i , as 

^P A D J , i = ^P SYSTEM  + U C T , i S^S  Y S T E M (-*) 

where  i -=  1 or  2 for  forward  or  backward  flow,  respectively.  Then,  for  a 
particular  value  of  ApF  there  could  be 

forward  flow  through  systems  with  non-zero  duct  resistance  if: 

(ApF  + ApADJ( x)  > 0 (6) 

Extending  the  ideas  presented  in  [1], 

when  Eq.  (6)  is  satisfied  and  there  is  forward  flow,  the  duct  flow 
resistance,  Rx , is  defined  as  the  value  which  would  lead  to  a mass  flow- 
rate  through  the  system  duct,  MDUCT  1 > 0,  and  corresponding  volume 
flow-rate,  VDUCT  x > 0,  which  satisfies: 

P-i  = (APp  + ApAD  j _ F ) 1 1 2/Mduct  > ! *=  (ApF  + ApAD  j ( j )1  / 2/VducT  . ! ^duct  , l 

(7) 
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Similarly,  there  could  be 

backward  flow  through  systems  with  non- zero  duct  resistance  if: 

(ApF  + ApADJ  2)  < 0 (8) 

and  when  Eq.  (8)  is  satisfied  and  there  is  backward  flow,  the  duct  flow 
resistance,  R2 , is  defined  as  the  value  which  would  lead  to  a mass  flow- 
rate  through  the  system  duct,  MDUCT  2 < 0,  and  corresponding  volume 
flow-rate,  VDUCT  2 < 0,  which  satisfies: 

R2  = '("APf  ■ APad j _ 2 ) 1 / 2/^DUCT , 2 = ‘(‘^Pf  ■ APaD j f 1 ) 1 ' 2/PdUCT , 2^DUCT , 2 

(9) 

Note  that  for  a particular  value  of  ApF  there  is  a possiblity  that  the 
conditions  of  Eqs . (6)  and  (8)  for  forward  and  backward  flow  can  be  satisfied 
simultaneously.  When  such  a circumstance  occurs  with  the  solution  value  of 
ApF , the  correct  one  of  the  two  possible  solutions  to  the  flow  problem  of  Eqs. 
(7)  and  (9),  would  require  other  considerations. 

If  the  duct  flow  resistance  is  negligible,  i.e.,  Rx  = R2  =0  then  ApF 
solutions  to 

ApF  = -ApADJ1  if  f(-ApADJ1)  > 0 

or  (10) 

aPf  " -aPadj,2  if  f ( “ aPAd j , 2 ) > 0 

and  the  corresponding  VQUCT  *=  VF  solutions  that  would  be  obtained  from  Eq.  (3) 
would  represent  valid  [VF , ApF ] solution  pairs  to  the  desired  fan/duct  flow 
problem. 


The  Problem  of  Determining  Possible  Solutions  for  the  System  Mass  Flow-Rate 

and  the  Corresponding  Apf  When  Duct  Resistance  is  Non-Zero 

For  a valid  fan/duct  system  flow  solution,  the  duct  volume  flow- rate  of  Eq. 
(7)  or  (9)  and  its  corresponding  mass  flow  rate  must  be  identical  to  the  fan 
volume  flow- rate  of  Eq.  (3)  and  its  corresponding  mass  flow- rate.  Equating 
these  volume  flow- rates  leads  to  the  following  result: 


FANRES 


3 


For  non-zero  duct  resistance,  ApF  solutions  to 


IaPf  + aPadj, 1 I1/2/(^duct, iRi)  " f(APF>  if  f(APF)  > 0 
or  (11) 


-|ApF  + ApADJ  2 |1/2/(pDUCt,2R2)  “ f(APp)  if  f(APp)  < 0 

and  the  corresponding  VDUCT  - V£  solutions  that  would  be  obtained  from 
Eq.  (3)  would  represent  valid  [VF , ApF ] solution  pairs  to  the  desired 
fan/duct  flow  problem. 

Define  the  VF  = 0 intercept  of  the  simulated  fan  cuve  of  Eq.  (3)  as  ApF  0 , 
i . e . , 


f (aPf  t o ) “ 0 

Then  it  can  be  shown  that: 


(12) 


A forward-flow  solution,  ApF  S0LN  , to  Eqs . 
only  if  ApADJ1  satisfies  -ApADJ1  < ApF  0 
bounded  by  -ApADJ1  < ApFS0LN  < ApF  0 . 


(10)  or  (11)  exists  if  and 
and  such  a solution  is 


(13) 


A backward- flow  solution,  ApF  S0LN  , to  Eqs. 
only  if  ApADJ2  satisfies  -ApADJ2  > ApF  „ , 
bounded  by  ApF  0 < ApFS0LN  < -ApADJ>2. 


(10)  or  (11)  exists  if  and 
and  such  a solution  is 


(14) 


It  also  follows  from  this  that 

double  solutions  to  Eqs.  (10)  or  (11)  require: 

APaDJ,1  " AP  A D J , 2 > 0 or  (PdUCT.I  ' U C T , 2 ^ A^  SYSTEM  > 0 

The  results  of  (13) -(15)  are  useful  in  the  determining  numerically  the 
aPf , s o l n solutions. 

If  Eqs.  (10)  or  (11)  do  not  lead  to  any  solutions,  then  the  solution  to  the 
duct/fan  flow  problem  is  one  of  zero  flow,  i.e., 

VDUCt  “ ^f  “ 0 and  ApF  - ApF  0 if  there  are  no  solutions  to 

Eqs.  (10)  or  (11) 


(16) 

There  are  four  possible  classes  of  solution  to  Eqs.  (10)  or  (11),  namely:  a 
single  forward- flow  solution  and  no  backward- flow  solution,  a single  backward- 
flow  solution  and  no  forward- flow  solution,  one  forward- flow  solution  and  one 
backward- flow  solutions,  and  no  solutions  at  all.  In  the  first  two  classes, 
the  unique  solution  obtained  is  taken  to  be  the  desired  solution  to  the 
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fan/duct  flow  problem.  When  two  solutions  are  obtained,  both  are  returned  as 
possible  solutions  to  the  fan/ duct  flow  problem. 

When  no  solutions  to  Eqs . (10)  or  (11)  are  obtained,  the  single  solution  of 
Eq.  (16)  is  taken  to  be  the  solution  to  the  fan/ duct  flow  problem. 


OUTPUT 


Msoln.n*  N = 1 or  2 

Solutions  for  the  mass  flow- rate  through  the  system. 

If  a flow  solution  is  from  from  room  1 to  room  2 
(normal  inlet  to  normal  outlet,  i.e.,  in  the  normal - 
fan-operation  flow  direction),  then  Mg0LN  N > 0.  If  a 
flow  is  from  room  2 to  room  1 (normal  outlet  to  normal 
inlet,  i.e.,  a condition  of  fan  backflow),  then 
Mg o l n N < 0 . Msoln  n “ 0 indicates  a zero  flow-rate 
solution.  If  there  is  only  one  solution,  i.e.,  Nsoln 
= 1,  then  it  is  Msqln.i-  tkS/s ] 


^SOLN  “ 1 °r  ^ 

Number  of  solutions. 


INPUT 


ACC 


Absolute  error  tolerance  for  the  fan/duct  system's 
computed  mass  flow  rate(s). 


VF  J,  J = 1 to  JTERM 

Values  of  the  volume  flow- rate  at  JTERM  points  used  to 
approximate  the  fan  operating  curve  (volume  flow- rate 
vs  cross-fan  pressure  difference).  [m3/s] 

JTERM 


Number  of  points  used  to  approximate  the  fan  operating 
curve . 
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Ri  [R2] 


Effective  resistance  of  the  duct  when  flow  is  in  the 
forward  [backward]  direction.  [Pa1/2kg_1s  ■= 

{kg/ (m-s2 ) } 1 1 2 kg' 1 s ] 


ApFJ1  J - 1 to  JTERM 

Values  of  the  cross -fan  pressure  difference  at  JTERM 
points  used  to  approximate  the  fan  operating  curve 
(volume  flow-rate  vs  cross-fan  pressure  difference). 
[Pa  *=  kg/(m-s2)] 

aPsystem 

Pressure  difference  across  the  system  end-points,  see 
Eq.  (2).  [Pa  = kg/(m-s2)] 

^SYSTEM 


Difference  in  elevation  of  the  system  endpoints,  see 
Eq.  (1).  [m] 


^DUCT , 1 t^DUCT , 2 ] 


Density  of  the  system  flow  when  flow  is  in  the  forward 
[backward]  direction,  [kg/m3 ] 


CALCULATIONS 

1.  Calculate  ApADJ  : and  ApADJ  2 from  Eq.  (5). 

2.  Calculate  the  forward- flow  solution  if  it  exists. 

3.  Calculate  the  backward- flow  solution  if  it  exists. 

4.  Return  to  the  calling  routine. 

a.  If  Mduct  y is  non-zero  and  MDUCT  2 - 0,  then  there  is  a unique 
forward  solution.  Set  Mg0LN  : “ mduct  i > nsoln  “ and  return  to 
the  calling  routine. 

b.  If  Mduct  1 - 0 and  MDUCT  2 non-zero,  then  there  is  a unique 
backward  solution.  Set  Ms0LN  x - ^duct  2»  ^solh  “ anc*  return 
to  the  calling  routine. 
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c.  If  MDUCT1  “ 0 and  MDUCT  2 “ 0.  then  there  is  a zero- flow 
solution.  Set  MS0LN  = 0,  Nsoln  = 1,  and  return  to  the  calling 
routine . 

d.  If  Mduct  i and  MDUCT  2 are  both  nonzero  then  there  are  two 
solutions,  one  for  forward  flow  and  one  for  backward  flow.  Set 

^duct.i  = ^duct,i>  ^soln,2  = ^duct,2>  nsoln  = 2>  and  return  to  the 
calling  routine. 


SUBROUTINES  USED 

INTERP  [3] 


REFERENCES 

[1]  Cooper,  L.Y.  and  Forney,  G.P.,  Section  3.4.5  of  Consolidated  Compartment 
Fire  Model  (CCFM)  Computer  Code  Application  CCFM. VENTS  - Part  I: 

Physical  Basis,  NISTIR  90-4342,  National  Institute  of  Standards  and 
Technology,  Gaithersburg  MD,  1990. 

[2]  Klote,  J.H.  and  Cooper,  L.Y. , Model  of  a Fan-Resistance  Ventilation 
System  and  Its  Application  to  Fire  Modeling,  NISTIR  89-4141,  National 
Institute  of  Standards  and  Technology,  Gaithersburg  MD,  September  1989. 

[3]  Forney,  G.P.,  Algorithm  INTERP  - Interpolates  A Table  of  Numbers, 
Appendix  B of  Consolidated  Compartment  Fire  Model  (CCFM)  Computer  Code 
Application  CCFM. VENTS  - Part  III:  Algorithms  and  Subroutines,  Cooper, 
L.Y. , and  Forney,  G.P.,  Eds.,  NISTIR  90-4344,  National  Institute  of 
Standards  and  Technology,  Gaithersburg  MD,  1990. 


SUBROUTINE  VARIABLES 

All  nomenclature  in  the  subroutine  is  identical  to  the  nomenclature  used 
above  except  for: 


MDUCT,l>  MD  U C T , 2 
^SOLN , 1 > ^SOLN , 2 
NS0LN 


aPf,  j 


XMFAN1 , XMFAN2  [kg/s] 

FANMAS (1) , FANMAS (2)  [kg/s] 

NSOLN 

FANVOL(J)  [m3/s] 

RESIS(l),  RESIS (2)  [Pa1/2kg_1s  = 
{kg/ (m*s2 ) )1/2kg"1s] 

FANPRS(J)  [Pa  - kg/(m-s2)] 
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A^  SYSTEM 

DYS  [m] 

aPad j , i > aPad J , 2 

DPFAJ  [Pa  - kg/(m-s2)] 

aPf,o 

DPFAN  [Pa  = kg/(m • s2 ) ] 

aPsystem 

DPS  [Pa  = kg/(m • s2 ) ] 

PdUCT,1»  PdUCT,2 

RHOSYS (1) , RHOSYS (2)  [kg/m3 

PREPARED  BY 

Leonard  Y.  Cooper,  Glenn  P. 

Forney,  and  John  H.  Klote 

July  1989 
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Figure  1.  The  fire -generated  environment  in  two  rooms  connected  by  a 
fan/ duct  system. 
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Figure  2. 


A sketch  of  a generic  fan- operating  curve  which  increases 
monotonically  in  volume -flow- rate  with  decreasing  cross-fan 
pressure  difference. 
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SUBROUTINE 


SUBROUTINE  FANRES ( 

I ACC , FANPRS , FANVOL , JTERM , DPS , RES IS , DYS , RHOSYS , 

0 FANMAS , NSOLN 

+ ) 

IMPLICIT  DOUBLE  PRECISION  (A-H.O-Z) 

C*BEG 

C***  FANRES  GENERIC  - FIND  THE  MASS  THAT  A FAN  WILL  DELIVER  UNDER 
C GIVEN  CONDITIONS  IN  THE  INLET  AND  OUTLET  ROOMS. 

C 

C***  SUBROUTINE  VARIABLES 


C 

C INPUT 

C 

C ACC 

C FANPRS 

C FANVOL 

C JTERM 

C DPS 

C RESIS 

C 

C DYS 

C RHOSYS 

C 

C OUTPUT 

C 

C FANMAS 

C NSOLN 

C*END 


REQUESTED  ACCURACY 

LIST  OF  JTERM  PRESSURES  [P] 

LIST  OF  JTERM  VOLUME  FLOWS  [M**3/S] 

NUMBER  OF  POINTS  IN  FAN  CURVE  DEFINED  BY  (FANPRS , FANVOL) 
PRESSURE  DIFFERENCE  ACROSS  THE  SYSTEM  END  POINTS  [P] 
RESITANCE  OF  FAN  SYSTEM  FOR  FORWARD  AND  REVERSE  FLOW 
CONDITIONS 

DIFFERENCE  IN  ELEVATION  OF  THE  SYSTEM  ENDPOINTS  [M] 
DENSITY  AT  INLET  AND  OUTLET  OF  FAN  [KG/M**3] 


MASS  DELIVERED  BY  FAN  [KG/S] 
NUMBER  OF  SOLUTIONS 


DIMENSION  FANPRS (*) .FANVOL (*) .RHOSYS (2) 

DIMENSION  FANMAS ( 2 ) , RESIS(2) 

EXTERNAL  FANFOR , FANREV 

RSISFW  - RESIS (1) 

RSISRV  •=  RESIS (2) 

C 

C***  GET  CANDIDATE  FORWARD  SOLUTION 
C 

C***  BRACKET  ROOT 
C 

DPMIN  -=  FANPRS (1) 

IF (JTERM . GE . 2 ) THEN 

DPMAX  - FANPRS (JTERM) 

ELSE 

DPMAX  = DPMIN  + 1.0D0 
END  IF 

IF(RSISFW.NE . 0 . 0D0)THEN 

CALL  BRACK (FANFOR, DPMIN, DPMAX, XA.XB, 

1 RSISFW , RHOSYS , DYS , DPS , FANPRS , FANVOL , JTERM) 

C 
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C***  REFINE  ESTIMATE  FOR  ROOT 
C 

CALL  RTSAF2 ( FANFOR , XA , XB , ACC , PFAN1 , 

1 RSISFW , RHOSYS , DYS , DPS , FANPRS , FANVOL , JTERM) 

C 

C***  CALCULATE  MASS  FLOW  THAT  CORRESPONDS  TO  ROOT 
C 

CALL  MRES ( PFAN1 , 1 , RSISFW , RHOSYS , DYS , DPS , XMFAN1 , XDUM) 
ENDIF 

IF(RSISRV . NE . 0 . 0D0)THEN 
C 

C***  GET  CANDIDATE  REVERSE  SOLUTION 

C 

C 

C***  BRACKET  ROOT 
C 

CALL  BRACK ( FANREV , DPMIN , DPMAX , XA , XB , 

1 RSISRV , RHOSYS , DYS , DPS , FANPRS , FANVOL, JTERM) 

C 

C***  REFINE  ESTIMATE  FOR  ROOT 
C 

CALL  RTSAF2 ( FANREV , XA , XB , ACC , PFAN2 , 

1 RSISRV , RHOSYS , DYS , DPS , FANPRS , FANVOL, JTERM) 

C 

C***  CALCULATE  MASS  FLOW  THAT  CORRESPONDS  TO  ROOT 
C 

CALL  MRES(PFAN2, -1, RSISRV, RHOSYS, DYS, DPS, XMFAN2, XDUM) 
ENDIF 

IF(RSISRV . EQ . 0 . 0D0 . OR . RSISFW . EQ . 0 . 0D0)THEN 
C 

C***  LOOK  FOR  SOLUTIONS  WHEN  RESISTANCE  OF  FAN  SYSTEM  IS  ZERO 
C 

G ■=  9.80D0 

DPFAJ  «=  DPS  + RHOSYS (1)*G*DYS 

CALL  INTERP ( FANPRS , FANVOL , JTERM , - DPFAJ , XMVOL , DXMVOL , 2 ) 
I F (XMVOL . GT . 0 . 0D0 ) THEN 

XMFAN1  - RHOSYS (l)*XMVOL 
ELSE 

XMFAN1  - 0.0D0 
ENDIF 

DPFAJ  - DPS  + RHOSYS ( 2 )*G*DYS 

CALL  INTERP ( FANPRS , FANVOL , JTERM , - DPFAJ , XMVOL , DXMVOL , 2 ) 
IF(XMVOL. LT . 0 . 0D0)THEN 

XMFAN2  - RHOSYS ( 2 )*XMVOL 
ELSE 

XMFAN2  -=  0.0D0 
ENDIF 
ENDIF 
C 

C***  DETERMINE  THE  VALID  SOLUTION(S) 

C 

IF (XMFAN1 . NE . 0 . 0D0 . AND . XMFAN2 . EQ . 0 . 0D0 ) THEN 
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c 

C***  A UNIQUE  SOLUTION  IN  THE  FORWARD  DIRECTION 
C 

FANMAS(l)  - XMFAN1 
NSOLN  -=  1 

ELSE  I F (XMFAN1 . EQ . 0 . 0D0 . AND . XMFAN2 . NE . 0 . 0D0 ) THEN 
C 

C***  A UNIQUE  SOLUTION  IN  THE  REVERSE  DIRECTION 
C 

FANMAS ( 1 ) -=  XMFAN2 
NSOLN  -=  1 

ELSE  IF (XMFAN1 . EQ . 0 . 0D0 . AND . XMFAN2 . EQ . 0 . 0D0 ) THEN 
C 

C***  A ZERO  FLOW  SOLUTION 
C 

FANMAS (1)  «=  0.0D0 
NSOLN  - 1 
ELSE 
C 

C***  TWO  SOLUTIONS 
C 

FANMAS (1)  = XMFAN1 
FANMAS (2)  = XMFAN2 
NSOLN  = 2 
ENDIF 
RETURN 
END 

SUBROUTINE  FANFOR ( DPFAN , XMNET , DXMNET , 

1 RES I S , RHO , DZ , DPS , FANPRS , FANVOL , NPTS ) 

IMPLICIT  DOUBLE  PRECISION  (A-H.O-Z) 

C*BEG 

C***  FANFOR  GENERIC  - CALCULATE  THE  NET  MASS  FLOW  IN  A FAN  SYSTEM 
C ASSUMING  THE  FLOW  IS  IN  A FORWARD  OR  NORMAL  DIRECTION. 

C THIS  ROUTINE  IS  CALLED  BY  A ZERO  FINDER  THAT  LOCATES  A 

C PRESSURE  DIFFERENCE  DPFAN,  THAT  CAUSES  THE  NET  MASS  FLOW 

C TO  BE  ZERO. 

C 

C***  SUBROUTINE  ARGUMENTS 
C 

C INPUT 


C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 


DPFAN 
RES  IS 
RHO 
DZ 
DPS 

FANPRS 

FANVOL 

NPTS 


OUTPUT 


PRESSURE  ACCROSS  FAN  [P] 

RESITANCE  OF  FAN  SYSTEM 

DENSITY  AT  INLET  AND  OUTLET  OF  FAN  [KG/M**3] 

DIFFERENCE  IN  ELEVATION  OF  THE  SYSTEM  ENDPOINTS  [M] 
PRESSURE  DIFFERENCE  ACROSS  THE  SYSTEM  END  POINTS  [P] 

LIST  OF  NPTS  PRESSURES  [P] 

LIST  OF  NPTS  VOLUME  FLOWS  [M**3/S] 

NUMBER  OF  POINTS  IN  FAN  CURVE  DEFINED  BY  (FANPRS , FANVOL) 
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C XMNET  NET  MASS  FLOW  [KG/(S*M**3) ] 
C DXMET  DERIVATIVE  OF  NET  FAN  FLOW 
C*END 


C 

C*** 

c 

c 

c*** 

c 

c 

c*** 

c 

c 


c 

c*** 

c 


C*BEG 

C*** 

C 

C 

c 

c 

c 

c*** 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


DIMENSION  RHO(2) , FANPRS (*) , FANVOL(*) 

CALCULATE  MASS  RESISTANCE  FLOW 
CALL  MRES (DPFAN , 1 , RESIS , RHO , DZ , DPS , XMRES , DXMRES ) 

CALCULATE  VOLUME  FLOW 

CALL  INTERP ( FANPRS , FANVOL , NPTS , DPFAN , XMVOL , DXMVOL , 2 ) 

CONVERT  VOLUME  FLOW  TO  MASS  FLOW  ASSUMING  FORWARD  DIRECTION 
(USE  DENSITY  IN  ROOM  1) 

XMFAN  = RHO ( 1 ) *XMVOL 
DXMFAN  - RHO ( 1 ) *DXMVOL 

COMPUTE  DIFFERENCE  BETWEEN  TWO  FLOWS 

XMNET  = XMFAN  - XMRES 
DXMNET  = DXMFAN  - DXMRES 
RETURN 
END 

SUBROUTINE  FANREV( DPFAN, XMNET, DXMNET, 

1 RESIS , RHO , DZ , DPS , FANPRS , FANVOL, NPTS ) 

IMPLICIT  DOUBLE  PRECISION  (A-H.O-Z) 

FANREV  GENERIC  - CALCULATE  THE  NET  MASS  FLOW  IN  A FAN  SYSTEM 
ASSUMING  THE  FLOW  IS  IN  A REVERSE  DIRECTION.  THIS 
ROUTINE  IS  CALLED  BY  A ZERO  FINDER  THAT  LOCATES  A 
PRESSURE  DIFFERENCE  DPFAN,  THAT  CAUSES  THE  NET  MASS 
FLOW  TO  BE  ZERO. 

SUBROUINE  VARIABLES 

INPUT 


DPFAN  PRESSURE  ACCROSS  FAN  [P] 

RESIS  RESITANCE  OF  FAN  SYSTEM 

RHO  DENSITY  AT  INLET  AND  OUTLET  OF  FAN  [KG/M**3] 

DZ  ELEVATION  OF  INLET  AND  OUTLET  OF  FAN  [M] 

DPS  PRESSURE  DIFFERENCE  ACROSS  THE  SYSTEM  END  POINTS  [P] 

FANPRS  LIST  OF  NPTS  PRESSURES  [P] 

FANVOL  LIST  OF  NPTS  VOLUME  FLOWS  [M**3/S] 

NPTS  NUMBER  OF  POINTS  IN  FAN  CURVE  DEFINED  BY  (FANPRS , FANVOL) 

OUTPUT 


XMNET  NET  MASS  FLOW  [M**3/S] 
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c 

C*END 


DXMET 


DERIVATIVE  OF  NET  FAN  FLOW 


DIMENSION  RHO(2) , FANPRS (*) , FANVOL(*) 

C 

C***  CALCULATE  MASS  RESISTANCE  FLOW 
C 

CALL  MRES(DPFAN, -1 ,RESIS ,RHO,DZ,DPS ,XMRES .DXMRES) 

C 

C***  CALCULATE  VOLUME  FLOW 
C 

CALL  INTERP ( FANPRS , FANVOL , NPTS , DPFAN , XMVOL , DXMVOL , 2 ) 

C 

C***  CONVERT  VOLUME  FLOW  TO  MASS  FLOW  ASSUMING  REVERSE  DIRECTION 
C (USE  DENSITY  IN  ROOM  2) 

C 

XMFAN  = RHO(2)*XMVOL 
DXMFAN  -=  RHO  ( 2 ) *DXMVOL 
C 

C***  COMPUTE  DIFFERENCE  BETWEEN  TWO  FLOWS 
C 

XMNET  = XMFAN  - XMRES 
DXMNET  = DXMFAN  - DXMRES 
RETURN 
END 

SUBROUTINE  MRES (DPFAN , IDIR , RESIS , RHO , DZ , DPS , XMRES , DXMRES ) 
IMPLICIT  DOUBLE  PRECISION  (A-H.O-Z) 

C*BEG 

C***  MRES  GENERIC  - THIS  ROUTINE  COMPUTES  THE  MASS  FLOW  IN  A FAN 
C SYSTEM  FOR  A GIVEN  RESITANCE  AND  PRESSURE 

C 

C***  SUBROUTINE  ARGUMENTS 
C 

C INPUT 

C 

C DPFAN 

C IDIR 

C 

C RESIS 

C RHO 

C DZ 

C DPS 

C 

C OUTPUT 

C 

C XMRES 

C DXMRES 

C*END 


PRESSURE  ACCROSS  FAN  [P] 

DIRECTION  OF  FLOW  IDIR>0  FOR  FORWARD  FLOW 
IDIR<0  FOR  REVERSE  FLOW 
RESISTANCE  OF  FAN  SYSTEM 

DENSITY  AT  INLET  AND  OUTLET  OF  FAN  [KG/M**3] 
DIFFERENCE  IN  ELEVATION  OF  THE  SYSTEM  ENDPOINTS  [M] 
PRESSURE  DIFFERENCE  ACROSS  THE  SYSTEM  END  POINTS  [P] 


MASS  RESISTANCE  FLOW  [KG/S] 
DERIVATIVE  OF  MASS  RESISTANCE  FLOW 


DIMENSION  RHO(2) , PROOM(2) 
PARAMETER  (G-9.800D0) 

I - 1 
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IF(IDIR. LT .0)1  - 2 
DPAJ  - DPS  + RHO(I)*G*DZ 
DPNET  - DPFAN  + DPAJ 
SQRNET  - SQRT(ABS (DPNET)) 

XMRES  = 0.0D0 
DXMRES  - O.ODO 

IF ( I . EQ . 1 . AND . DPNET . GT . 0 . ODO ) THEN 
XMRES  - SQRNET/RESIS 
DXMRES  = . 5D0/(SQRNET*RESIS) 

ELSE  IF ( I . EQ . 2 . AND . DPNET . LT . 0 . ODO ) THEN 
XMRES  -=  -SQRNET/RESIS 
DXMRES  - . 5D0/ ( SQRNET*RES I S ) 

ENDIF 

RETURN 

END 

SUBROUTINE  BRACK ( FUNCD , A , B , XA , XB , 

1 RES I S , RHO , DZ , DP , FANPRS , FANVOL , NPTS ) 

IMPLICIT  DOUBLE  PRECISION  (A-H.O-Z) 

C*BEG 

C 


C*** 

c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 


BRACK  SPECIFIC  - THIS  ROUTINE  BRACKETS  THE  ROOT  OF  THE 

FUNCTION,  FUNCD  SO  THAT  THE  ROOT  LIES  IN  THE  INTERVAL 
[XA,XB] . THIS  SUBROUTINE  IS  BASED  ON  A SUBROUTINE 
NAMED  ZBRAC  FOUND  IN: 

NUMERICAL  RECIPES,  THE  ART  OF  SCIENTIFIC  COMPUTING,  PRESS ; 
WILLIAM  H.,  FLANNERY,  BRAIN  P.,  TEUKOLSKY,  SAUL  AND 
VETTERLING,  WILLIAM  T. , CAMBRIDGE,  PP.  245,  1986 

SUBROUTINE  ARGUMENTS 


INPUT 

FUNCD 

A 

B 

RES  IS 
RHO 
DZ 
DP 

FANPRS 

FANVOL 

NPTS 

OUTPUT 

XA 

XB 


FUNCTION  TO  BE  BRACKETED 
CANDIDATE  LEFT  BRACKET 
CANDIDATE  RIGHT  BRACKET 
RESITANCE  OF  FAN  SYSTEM 

DENSITY  AT  INLET  AND  OUTLET  OF  FAN  [KG/M**3] 

DIFFERENCE  IN  ELEVATION  OF  THE  SYSTEM  ENDPOINTS  [M] 
PRESSURE  DIFFERENCE  ACROSS  THE  SYSTEM  END  POINTS  [P] 

LIST  OF  NPTS  PRESSURES  [P] 

LIST  OF  NPTS  VOLUME  FLOWS  [M**3/S] 

NUMBER  OF  POINTS  IN  FAN  CURVE  DEFINED  BY  (FANPRS , FANVOL) 


ACTUAL  LEFT  BRACKET 
ACTUAL  RIGHT  BRACKET 


C*END 


DIMENSION  RHO(*) , FANPRS (*) , FANVOL(*) 
PARAMETER  (MAXIT-50) 
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AL  = A 
AU  = A 
BL  ■=  B 
BU  = B 

CALL  FUNCD (AL , FAL , DFAL , RES IS , RHO , DZ , DP , FANPRS , FANVOL, NPTS ) 
CALL  FUNCD ( BL , FBL , DFBL , RES I S , RHO , DZ , DP , FANPRS , FANVOL , NPTS ) 

FBU  = FBL 
FAU  = FAL 

DO  10  I = 1,  MAXIT 

IF ( FAL* FBL . LE . 0 . 0D0 ) THEN 
XA  = AL 
XB  = BL 
RETURN 
END  IF 
AU  = BU 
FAU  = FBU 

BU  = BU  + (BU-A)*2.0D0 

CALL  FUNCD ( BU , FBU , DFBU , RES I S , RHO , DZ , DP , FANPRS , FANVOL , NPTS ) 
I F ( FAU* FBU . LE . 0 . 0D0 ) THEN 
XA  = AU 
XB  - BU 
RETURN 
END  IF 
BL  = AL 
FBL  = FAL 

AL  = AL  - (B-AL)*2 . 0D0 

CALL  FUNCD ( AL , FAL , DFAL , RES I S , RHO , DZ , DP , FANPRS , FANVOL , NPTS ) 
10  CONTINUE 

CALL  MSGPRT ( ' ROOT  NOT  BRACKETED  IN  BRAKIT',3) 

RETURN 

END 

SUBROUTINE  RTSAF2 ( FUNCD , XI , X2 , XACC , RTSAFE , 

1 RES I S , RHO , DZ , DP , FANPRS , FANVOL , NPTS ) 

IMPLICIT  DOUBLE  PRECISION  (A-H.O-Z) 

C*BEG 

C 

C***  RTSAF2  SPECIFIC  - THIS  ROUTINE  FINDS  A VALUE  OF  RTSAFE  SUCH 
C THAT  FUNCD  (RTSAFE)  *=  0.  IT  USES  A MODIFICATION  OF 

C NEWTON'S  METHOD.  THIS  SUBROUTINE  IS  BASED  ON 

C A SUBROUTINE  NAMED  RTSAF  FOUND  IN: 


C 

C 

c 

c 

c 


NUMERICAL  RECIPES,  THE  ART  OF  SCIENTIFIC  COMPUTING,  PRESS, 
WILLIAM  H.,  FLANNERY,  BRAIN  P.,  TEUKOLSKY,  SAUL  AND 
VETTERLING,  WILLIAM  T.,  CAMBRIDGE,  P.  258,  1986 


C***  SU] 
C 

C INPUT 


SUBROUTINE  ARGUMENTS 


C 

C FUNCD 
C XI 
C X2 


FUNCTION  WHOSE  ZERO  IS  DESIRED 
LEFT  BRACKET 
RIGHT  BRACKET 
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c 

XACC 

ACCURRACY  REQUEST 

c 

RES  IS 

RESITANCE  OF  FAN  SYSTEM 

c 

RHO 

DENSITY  AT  INLET  AND  OUTLET  OF  FAN  [KG/M**3] 

c 

DZ 

DIFFERENCE  IN  ELEVATION  OF  THE  SYSTEM 

ENDPOINTS 

[M] 

c 

DPS 

PRESSURE  DIFFERENCE  ACROSS  THE  SYSTEM 

END  POINTS 

[P] 

c 

FANPRS 

LIST  OF  NPTS  PRESSURES  [P] 

c 

FANVOL 

LIST  OF  NPTS  VOLUME  FLOWS  [M**3/S] 

c 

NPTS 

NUMBER  OF  POINTS  IN  FAN  CURVE  DEFINED 

BY  (FANPRS .FANVOL) 

c 

C OUTPUT 

C 

C 

C RTSAFE  ZERO  OF  THE  FUNCTION  FUNCD  WITHIN  A TOLERANCE  OF  XACC 
C*END 


DIMENSION  RHO(2) ,FANPRS(*) ,FANVOL(*) 

PARAMETER  (MAXIT-100) 

CALL  FUNCD (XI , FL , DF , RES IS , RHO , DZ , DP , FANPRS , FANVOL , NPTS ) 
IF(FL. EQ. 0 . 0D0)THEN 
RTSAFE  = XI 
GO  TO  999 
END  IF 

CALL  FUNCD (X2 , FH , DF , RES IS , RHO , DZ , DP , FANPRS , FANVOL , NPTS ) 
IF(FH . EQ . 0 . 0D0)THEN 
RTSAFE  - X2 
GO  TO  999 
END  IF 

IF(FL*FH.GT. 0 . )THEN 

CALL  MSGPRT( ' ROOT  MUST  BE  BRACKETED  IN  RTSAF2',3) 
RETURN 
END  IF 

IF ( FL . LT . 0 . ) THEN 
XL— XI 
XH-X2 
ELSE 
XH-X1 
XI^X2 
SWAP-FL 
FL-FH 
FH-SWAP 
END  IF 

RTSAFE-. 5* (X1+X2) 

DXOLD— ABS (X2 -XI ) 

DX-DXOLD 

CALL  FUNCD (RTSAFE , F , DF , RES IS , RHO , DZ , DP , FANPRS , FANVOL , NPTS ) 
DO  11  J-l.MAXIT 

IF( ( (RTSAFE -XH)*DF-F)*( (RTSAFE -XL) *DF-F) .GE.O. 

* .OR.  ABS (2 . *F) . GT . ABS (DXOLD*DF)  ) THEN 

DXOLD-DX 
DX-0 . 5*(XH-XL) 

RTSAFE-XL+DX 
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IF(XL. EQ.RTSAFE)GO  TO  999 
ELSE 

DXOLD-DX 
DX-F/DF 
TEMP-RTSAFE 
RTSAFE=RTSAFE - DX 
IF (TEMP. EQ.RTSAFE)GO  TO  999 
END  IF 

IF(ABS(DX) .LT.XACC)GO  TO  999 

CALL  FUNCD (RTSAFE , F , DF , RES IS , RHO , DZ , DP , FANPRS , FANVOL , NPTS ) 
IF(F. LT . 0 . ) THEN 
XL-RTSAFE 
FL=F 
ELSE 

XH=RTSAFE 
FH=F 
END  IF 

11  CONTINUE 

CALL  MSGPRT(  ' NON -CONVERGENCE  IN  RTSAF2\3) 

RETURN 
999  CONTINUE 
RETURN 
END 
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FL0G02 


Deposition  of  Mass,  Enthalpy,  Oxygen,  and  Other  Product-of- 
Combustion  Flows  Passing  Between  Two  Rooms  Through  a Vertical, 
Constant -Width  Vent  or  Through  a Forced  Ventilation  System 


DESCRIPTION 


Consider  an  instant  of  time  during  the  simulation  of  a multi -room  compartment 
fire  environment.  The  flow  through  a constant-width  natural  vent  located  in  a 
wall  segment  common  to  any  two  rooms  of  the  compartment  can  be  characterized 
generally  as  having  NSLAB  < 6,  uniform-property,  locally-near-horizontal , 
adjacent  slabs  of  flow.  This  is  depicted  in  Figure  1.  The  elevations  of  the 
top  and  bottom  of  each  of  these  flow  slabs  is  such  that  it  is  deposited 
entirely  into  only  one  of  the  layers  of  the  receiving  room.  Also,  the  source 
of  each  flow  slab  is  entirely  from  one  of  the  two  layers  of  the  source  room. 

The  flow  through  a forced- ventilation  system  connecting  two  fire  - generated 
room  environments  can  often  be  characterized  as  being  fully-mixed  and  uniform- 
property.  In  the  duct  and  near  the  inlet  vent  this  flow  would  be  formed  by 
one  or  two  flow  slabs  entering  the  vent  from  the  source  room,  one  flowing  from 
the  upper  layer  and  one  flowing  from  the  lower  layer.  This  is  depicted  in 
Figure  2 for  the  case  of  two  flow  slabs.  The  duct  flow  would  be  ejected  at 
the  outlet  vent  into  the  recieving  room  as  one  or  two  identical -property  flow 
slabs , one  flow  slab  entering  the  upper  layer  and  one  entering  the  lower 
layer . 

In  the  case  of  the  unforced  natural  vent,  each  of  the  flow  slabs  leaves  one 
room  and  enters  the  immediately  adjacent  room  at  the  same  elevation  interval. 
In  the  case  of  the  forced-ventilation  system,  a flow  slab  generally  leaves  one 
room  and  enters  the  other  at  two  completely  different  elevations.  Indeed,  the 
two  rooms  joined  by  the  ventilation  system  will  generally  be  displaced 
arbitrarily  from  one  another  both  horizontally  and  vertically. 

It  is  assumed  that  a complete  description  of  the  flow  slabs  for  the  time  and 
natural  vent  or  forced-ventilation  system  of  interest  has  been  obtained  and  is 
available.  The  present  algorithm  determines  the  contributions  of  the  flow 
slabs  to  the  instantaneous  rates  of  flow  of  mass,  enthalpy,  oxygen  and  other 
products  of  combustion  into  the  upper  and  lower  layer  of  each  of  the  two  rooms 
being  affected  [1].  When  this  algorithm  is  used  for  forced  ventilation  flows 
a separate  application  of  the  algorithm  would  be  required  for  each  inlet  and 
outlet  vent.  In  the  case  of  an  inlet  vent,  the  modifications  to  the  source - 
room  layers  would  be  calculated  by  the  algorithm  and  implemented,  while  the 
calculated  effects  to  the  receiving  room  would  be  ignored,  i.e.,  the  receiving 
room  would  be  treated  as  a "dummy"  room.  In  the  case  of  an  outlet  vent,  the 
modifications  to  the  receiving -room  layers  would  be  calculated  and  implemented 
while  those  to  the  source  room  would  be  ignored. 

The  input  descriptions  of  the  flow  slabs  is  assumed  here  to  include  their 
number,  NSLAB , and  for  each  of  these:  the  direction  of  its  flow,  i.e.,  into 
room  1 or  room  2,  as _ identified  by  IDIR  N,  total  mass  flow  rate,  MsLAB,n> 
enthalpy  flow  rate,  QSLAB  N,  mass  flow  rate  of  oxygen,  P02  SLAB  N,  flow  rates 
of  NpR0D  - 1 other  products  of  combustion  being  tracked  by’ the  simulation, 
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PKslab,n>  K - 2 to  NpR0D , effective  heights  of  the  flow  slabs  in  rooms  1 and 
2 above  a datum  elevation,  ySLAB  x N and  ySLAB  2 n > anc^  absolute 
temperature,  TSLAB  N,  of  the  slab  upon  penetrating  the  receiving  room,  after 
deceleration  to  relatively  small  velocity,  but  before  significant  mixing. 

Also  assumed  to  be  known  in  each  of  the  two  rooms  is  the  absolute  temperature 
of  the  upper  and  lower  layers  (if  they  are  of  nonzero  thickness),  Ty  z and 
Tl  j , and  the  elevation  relative  to  the  datum  elevation  of  the  ceiling,  layer 
interface  and  floor,  yCEIL>I,  yLAYER>I,  and  yFLOOR.i-  1 - 1 or  2 , 
respectively.  When  one  is  dealing  with  slab  flows  through  an  unforced 
vertical  vent  in  a common  wall  between  adjacent  rooms  the  input  data  should 
satisfy  ySLAB>1>H  - ySLAB.2,N  for  a11  Nslab  slabs.  As  discussed  above,  this 
is  not  generally  true  for  the  flow  through  a forced-ventilation  system  where 
the  system  endpoint  elevations  are  generally  different. 

FL0G02  determines  the  source  and  destination  layer (s)  of  all  the  slab  flows  of 
a vent  and  the  net  contributions  of  these  slab  flows  to  the  rates  of  increase 
of  mass,  My  z and  ML  r,  enthalpy,  Qy  z and  QL  and  products  of  combustion, 

PK  u i and  PK  L j , of  the  upper  and  lower  layers  of  the  two  rooms.  The  en- 
thalpy-flow-rate calculation  includes  heat  transfer  exchanges  between  the 
incoming  enthalpy  flows  and  the  bounding  surfaces  of  the  receiving  rooms.  The 
algorithm  for  mass  apportionment  of  a flow  slab  entering  the  recieving  room  is 
based  on  an  extension  to  a simple  scheme,  proposed  on  pages  24-25  of  [2],  for 
depositing  into  or  apportioning  between  the  upper  and  lower  layers  of  the 
receiving  room,  the  flow  of  any  particular  slab.  In  determining  the  net 
enthalpy  flow  rate  to  a receiving  layer  associated  with  a given  mass  flow  rate 
contribution,  an  algorithm  based  on  an  extension  of  ideas  proposed  in  [3]  and 
used  in  [4]  for  heat  transfer  losses  to  the  walls  of  the  room  is  adopted. 

The  rule  of  apportioning  the  entering  slab  flow  has  the  following  specific 
elements : 

1.  If  a)  the  temperature  of  the  slab  is  greater  than  the  temperature 
of  the  lower  layer  in  the  receiving  room,  but  by  an  amount  less 
than  a small  user- specified  value  ATEps , and  b)  the  upper  layer  of 
the  receiving  room  is  of  zero  thickness,  i.e.,  the  layer  interface 
is  at  the  ceiling  of  the  receiving  room,  then  c)  the  entire  slab 
is  deposited  in  the  lower  layer  and  no  new  upper  layer  is  allowed 
to  grow  by  virtue  of  mass  additions  from  the  slab  flow. 

2.  If  a)  the  temperature  of  the  slab  is  less  than  the  temperature  of 
the  upper  layer  in  the  receiving  room,  but  by  an  amount  less  than 
a small  user- specif ied  value  ATEps , and  b)  the  lower  layer  of  the 
receiving  room  is  of  zero  thickness,  i.e.,  the  layer  interface  is 
at  the  floor  of  the  receiving  room,  then  c)  the  entire  slab  is 
deposited  in  the  upper  layer  and  no  new  lower  layer  is  allowed  to 
grow  by  virtue  of  mass  additions  from  the  slab  flow. 

3.  If  neither  of  the  above  two  elements  are  relevant,  and  the 
temperature  of  a slab  is  greater  than  the  temperature  of  the 
receiving  room  upper  layer,  then  the  entire  slab  flow  is  deposited 
in  the  upper  layer.  If  it  is  less  then  that  of  the  lower  layer, 
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then  the  entire  slab  flow, is  deposited  in  the  lower  layer.  If  the 
slab  temperature  is  between  the  temperature  of  the  upper  and  lower 
layer  then  the  slab  flow  is  divided  between  the  two  layers  in 
direct  proportion  to  the  temperature  differences. 

The  rule  of  determining  losses  of  slab  enthalpy  flow  by  means  of  heat  transfer 
to  the  upper  or  lower  surfaces  of  the  receiving  room  has  the  following 
specific  elements: 

1.  The  "excess"  enthalpy  flow  rate  of  a slab,  calculated  relative  to 
the  temperature  of  the  layer  which  receives  it  and  the  temperature 
of  the  slab  itself,  is  determined.  This  is  apportioned  along  with 
the  "base"  enthalpy  flow  rate,  calculated  at  the  temperature  of 
the  receiving  layer,  to  the  upper  layer  and/or  to  the  lower  layer, 
in  a manner  consistent  with  the  above  rule  for  apportioning  the 
mass  flow 

2.  A user-specified  fraction,  Av , of  the  upper-  and  lower-layer- 
portions  of  the  "excess"  enthalpy  flow  rate  is  assumed  to  be 
transferred  to  the  upper  and  lower  layer  surfaces  of  the  receiving 
room.  The  remaining  enthalpy  flow  rates  are  added  to  the  upper 
and  lower  layer. 

Although  the  following  calculations  are  not  actually  implemented  in  the 
present  subroutine,  it  is  noted  that  there  are  formulations  of  the  basic  model 
equations  where  it  would  be  useful  to  determine  explicitly  the  individual  and 
net  contributions  of  the  flow  slabs  to  the  group  of  terms 

Vi  " (CpVjiy  x - QU(I) 

Vi  = <CPTL(IML>I  - Ql>i) 

in  both  rooms  I.  Reference  here  is  to  formulations  which  use  layer  densities 
or  temperatures  as  independent  variables  in  the  differential  equations 
governing  the  room  environments . Under  such  circumstances  the  above  groups 
not  only  appear  explicitly  as  source  terms  in  the  generic  equation  set,  but 
they  often  have  the  property  of  being  zero  identically,  even  when  the  separate 
My  j and  Q0  z or  ML  : and  QL  x components  of  such  slab  flows  are  non-zero. 

For  example,  in  its  source  room  the  contribution  to  these  groups  of  each  flow 
slab  will  be  zero  identically.  When  required,  the  calculation  of  the  Gy  x and 
Gl  j groups  can  be  carried  out  in  a manner  as  to  insure  that  this  fact  is 
properly  taken  account  of,  i.e.,  to  insure  that  unneccessary  introduction  of 
roundoff  error  in  the  calculation  of  these  groups  is  eliminated. 


OUTPUT 


Vi  [Gl  x];  I - 1 or  2 

(NOTE:  THE  PRESENT  SUBROUTINE  DOES  NOT  PRODUCE  THIS  OUTPUT) 

Total  contribution  of  Gy  x - (CpTy  jMy  z - Qy _ z ) [GL  x 
- (CpTL  jMl  j - Ql  j)]  to  the  upper  [lower]  layer  of 
room  I due  to  all  N$LAB  slabs  of  the  vent  flow  [W] . 
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MU(i  [Mli];  I - 1 or  2 


Total  rate  of  mass  flow  to  the  upper  [lower]  layer  of 
room  I due  to  all  NSLAB  slabs  of  the  vent  flow  [kg/s]. 


^0  2 , U , I [ ^0  2 , L , I 1 ' ^ ^ or  ^ 

Total  rate  of  oxygen  mass  flow  to  the  upper  [lower] 
layer  of  room  I due  to  all  N<.,  4B  slabs  of  the  vent 
flow  [ (kg  of  02)/s] . 

^k.u.i  K " 2 to  nprod  * I - 1 or  2 

Total  rate  of  flow  of  product  of  combustion  K to  the 
upper  [lower]  layer  of  room  I due  to  all  NSLAB  slabs 
of  the  vent  flow  [(unit  of  product  K)/s] . 

Qu.i  IQl.i1;  I - 1 or  2 

Total  rate  of  enthalpy  flow,  based  on  absolute 
temperature,  to  the  upper  [lower]  layer  of  room  I due 
to  all  Nslab  slabs  of  the  vent  flow  [W] . 


Qs  U , I [Qsl-.i1;  I “ 1 or  2 

Total  rate  of  heat  transfer  to  the  upper  [lower] 
bounding  surfaces  of  room  I due  to  all  NSLAB  slabs  of 
the  vent  flow  [W] . 


INPUT 


Specific  heat  at  constant  pressure  of  the  vent  flow 
[W-s/(kg-K)]  (suggest  103  W-s/(kg-K)  for  air  as 
default) . 


XD I R , N I N - 1 to  NSLAB 

A measure  of  the  direction  of  the  room  1-to-room  2 
flow  of  slab  N.  For  flow  from  room  1 to  room  2 or 
from  room  2 to  room  1,  the  value  of  IDIR  N is  1 or  -1, 
respectively.  IDIR  N is  0 if  the  room- to -room  flow  is 
identically  zero. 


FL0G02  - 4 


^SLAB,N>  N “=  1 tO  Nslab 

Total  mass  flow  rate  in  slab  N [kg/s] . 

nprod 

number  of  products  of  combustion,  being  tracked  in  the 
simulation. 


number  of  slabs  of  flow,  < 6,  between  the  bottom  and 
top  of  the  vent. 

*02, SLAB, N1  N = 1 tO  Nslab 

Mass  flow  rate  of  oxygen  in  slab  N [(kg  of  02)/s]. 

**k  , s l a b , n ’ K *=  2 to  NpR0D  ; N » 1 to  Nslab 

Flow  rate  of  product  of  combustion  K in  slab  N [ (unit 
of  product  K)/s] . 

Qslab,n>  N ■=  1 to  Nslab 

Total  enthalpy  flow  rate  in  slab  N based  on  absolute 
temperature  [W] . 


1 to  No 


Absolute  temperature  of  the  flow  in  slab  N as  it 
enters  the  receiving  room  [K] . 


Absolute  temperature  of  the  upper  [lower]  layer  in 
room  I [K]  . 
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ySLAB,l,N  [ySLAB,2,lJ>  N“  1 t0  NSLAB 

Effective  elevation  of  slab  N in  room  1 [room  2]  above 
the  datum  elevation.  In  cases  of  slab  flows  through 
natural  vents  between  adjacent  rooms  the  input  data 
should  satisfy  ySLAB>1>N  - ysL  ab  , 2 , n ■ This  is  not 
generally  true  for  the  flow  through  a ventilation 
system,  [m] 

yFLooR.i  tycEiL.i]  Flayer, i^»  I - 1 or  2 

Elevation  of  the  floor  [ceiling]  {layer  interface)  in 
room  I above  the  datum  elevation  [m] . 


ATp  p g 

Minimum  temperature  difference,  > 0,  between  TSLAB  N 
and  TL  j [Tjj  j ] required  to  initiate  growth  of  an 
upper  [lower]  layer  (in  the  recieving  room  I)  from  a 
zero  thickness  [K] . 


Av 


Fraction  of  the  source  slab's  rate -of -enthalpy  inflow, 
computed  relative  to  the  temperature  in  the  receiving 
room  at  the  effective  elevation  of  the  slab,  which  is 
assumed  to  be  lost  to  the  bounding  walls  of  the 
receiving  room. 


CALCULATIONS 

1.  Set  equal  to  zero  all  HV1,  hLl,  P02U(I,  P02,l(i«  , u . 1 « 

pk , l , i » Qu , i > Ql , i * Qsu , i » Qsl , i * K “ 2 to  Nprod ,1-1  and  2 . 

2.  Consider  each  of  the  NSLAB  flow  slabs  and 

a.  Determine  the  destination  (upper,  lower  or  both  layers  of  room  1 

or  2)  of  any  particular  slab,  slab  N.  Let  room  I be  the  receiving 
room.  Using  the  temperature  and  the  mass,  enthalpy,  oxygen,  and 
product  flow  rates  of  slab  N,  adjust  Gy  j,  GL  j,  My  x,  ML  j , 

P0  2 , U , I » P0  2 , L , I ’ Qu  , I > Ql  , I > Qs  U . I > Qs  L , I * and  PK  , U , I ' PK  , L , I ’ K 
- 2 to  NpR0D  according  to  the  flow  diagram  of  Figure  3. 
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b. 


Determine  the  origin  (upper  or  lower  layer  of  room  1 or  2)  of  slab 
N.  Let  room  J be  the  room  of  origin.  Then  according  to  the 
temperature  and  the  mass,  enthalpy,  oxygen  and  other  product  flow 
rates  of  slab  N,  adjust  appropriately  M0  ]t  ML  j,  P02  w Jt 
P0  2,l,J’  , J > Ql  , j > and  pk  , u , j > V l , J 1 K = ^ to  NpR0D  . Do  not 

adjust  Gy  j , Gl  j,  Qsu  j,  Qsl  j.  There  are  two  possible  cases. 


Case  Ila: 


ysLAB,j,N  > Flayer, j 


Slab  flow  comes  from  upper  layer  of  room  of  origin. 
Subtract  M$LAB>N  from  Mu  j , P02 , slab , n . from  ^02,u,j> 
QsLAB.N  from  Qu.J’  and.PK,SLAB.N  f rom . PK , U . J - K = 2 to 
NpROD-  Do  not  adjust  Mlj,  P02,l,j>  Ql  . j > or  the 


Case  lib: 


ysLAB,j,N  < Flayer , j 


Slab  flow  comes  from  lower  layer  of  room  of  origin. 
Subtract  MSLAB>N  from  ML (J , P02 , slab , n . frora  po2,l,j, 
QsLAB.N  from  Ql.J.  and.PK,SLAB,N  f r°m . PK . U , J - K = 2 to 
Nprod-  Do  not  adjust  Mu  j,  P02jU(J,  Qu , j , or  the 

PK  , U . J ' S • 
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SUBROUTINE  VARIABLES 

CP 

PU>  PL 
PDIR  . N 

. I ■ . 2 


CP  [W- s/ (kg-K) ] 

FU,  FL 
DIRS12 (N) 

UFLW2 (1,3,1) , UFLW2 (2 ,3,1)  [W] 
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M 

USLAB,N 

- 

XMSLAB(N)  [kg/s] 

M M 

A1U  , 1 > 1AU  , 2 

- 

UFLW2 (1,1,2) , UFLW2 (2 ,1,2)  [kg/s] 

Ml.l  Ml>2 

- 

UFLW2 (1,1,1),  UFLW2 (2,1,1)  [kg/s] 

N 

i’  P R 0 D 

- 

NPROD 

nslab 

- 

NSLAB 

p 

r0  2 , SLAB  , N 

- 

PSLAB(1,N)  [(kg  of  02)/s] 

P0  2 , U , 1 > P0  2 , U , 2 

- 

UFLW2 (1 ,3,2),  UFLW2 (2,3,2)  [(kg  of  02 

)/s 

p P 

r0  2 , L , 1 ’ r0  2 , L , 2 

- 

UFLW2 (1,3,1) , UFLW2 (2,3,1)  [(kg  of  02 

)/s 

, S L A B , N 

- 

PSLAB(K,N)  [(unit  of  product  K)/s] 

PK,U,1’  PK , U , 2 

- 

UFLW2 ( 1 , 2+K , 2 ) , UFLW2(2 , 2+K, 2)  [(unit 
product  K)/s] 

of 

PK  , L , 1 ’ PK  , L , 2 

- 

UFLW2 ( 1 , 2+K , 1 ) , UFLW2 ( 2 , 2+K , 1 ) [(unit 
product  K)/s] 

of 

Qs  L A B , N 

- 

QSLAB(N)  [W] 

Qu,l>  Qu,2 

- 

UFLW2 (1 , 2,2),  UFLW2 (2,2,2)  [W] 

Ql,1-  Ql,2 

- 

UFLW2 (1,2,1) , UFLW2 (2,2,1)  [W] 

Qsu  , 1 > Qs  U , 2 

- 

QSUVNT(l) , QSUVNT(2)  [W] 

Qs  L , 1 ’ Qs  L , 2 

- 

QSLVNT(l) , QSLVNT(2)  [W] 

T 

XSLAB  , N 

- 

TSLAB(N)  [K] 

Tu-i.  TU  2 

- 

TU(1),  TU(2)  [K] 

TL.l.  TL,2 

- 

TL(1) , TL(2)  [K] 

yCEIL , 1 > yCEIL , 2 

- 

YCEIL(l),  YCEIL(2)  [m] 

yFLOOR , 1 > yFLOOR , 2 

- 

YFLOR(l),  YFLOR ( 2 ) [m] 

yLAYER , 1 * yL AYER , 2 

- 

YLAY(l),  YLAY(2)  [m] 

ySLAB  , 1 , N > ySLAB  , 2 , N 

- 

YSLAB1 (N) , YSLAB2 (N)  [m] 

ateps 

- 

TEPS  [K] 

Av 

- 

RLAM 
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Figure  1 . 


Possible  distributions  of  the  cross-vent  pressure -difference 
distribution,  Apx  2 , for  a single  natural  vent  and  the  flow 
directions  of  its  one-to-six  uniform-property  flow  slabs. 
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LAYER,! 


SYSTEM,! 


y 

FLOOR,! 

i L 

sssssssssssss* 


Figure  2 Two  flow  slabs  entering  the  inlet  vent  of  the  duct  of  a forced- 
ventilation  system. 
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Figure  3. 
(Page  1 
of  2) 


Flow  Diagram  for  FL0G02  Adjustments  to  Gy  j , 

P0  2 , U , I * P0  2 , L , I » PK  . U , I * PK  . L , I » , I * , I 

Receiving  Room  I Due  to  a Vent  Flow  Slab 


. i * , I > . I • 

. Qso.1 » Qsl.1  in  the 
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Set  Fl  «=  1 - Fy  and  add: 

FU  *^SLAB  ,N  t0  , I 

FU ' F02 , SLAB , H t0  F02,U,I 

fu'fk,slab.n  t0  K , U , I * K“l.  NPROD 

FU  "OsLAB  , H ’ ■^■r/'^SLAB  , H ^ t0  Qu  , I 

FU  *Cp  ’ ^SLAB . R ^ [T„ , x “TSLAB , H 1 

+*v [tslab . h -Tr 1 ) 

to  G0  I 

‘FU  "QsLAB  . N ‘ [F_Tr/^-SLAB  . H 1 t0  QsU  , I 
^'MsLAB.H  t0  ^L.I 
FL * F02 , SLAB . N t0  F02 , L , I 
FL*FK,SLAB,R  to  FK,L,I»  K”F  to  NPROD 
FL  ‘QsLAB , H ' tl'^V tl- "Tr/TsLAB , N n to  Ql.I 
FL  *Cp  ’^SLAB . H t tTL . I “TSLAB , H 1 

+*V lTSLAB . H 'TR 1 ) 

tO  GL  j 

*FL "QsLAB , N * [T. "Tr/Tslab ,H ] to  Qsl.i 


Figure  3 
(Page  2 
of  2) 


Flow  Diagram  for  FL0G02  Adjustments  to  Gy  j , 

F02,D,I»  F0 2 , L , I * FK,U,I»  FK,L,I»  Qu , I » Ql  , I 

Receiving  Room  I Due  to  a Vent  Flow  Slab 


, I » „ I » ^L  . i > 

• Qsu.i.  Qsl.i  in  the 
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non 


SUBROUTINE 


SUBROUTINE  FL0G02( 

I DIRS 12 , YSLAB1 , YSLAB2 , XMSLAB , TSLAB, NSLAB , TU , TL , YFLOR , 

I YCEIL , YLAY , QSLAB , RLAM , 

I PSLAB , MXPRD , NPROD , MXSLAB , DTEPS , 

0 QSLVNT , QSUVNT , UFLW2 ) 

IMPLICIT  DOUBLE  PRECISION  (A-H.O-Z) 

C*BEG 

C***  FL0G02  GENERIC  - DEPOSITION  OF  MASS,  ENTHALPY,  OXYGEN,  AND  OTHER 
C PRODUCT - OF- COMBUSTION  FLOWS  PASSING  BETWEEN  TWO  ROOMS 

C THROUGH  A VERTICAL,  CONSTANT-WIDTH  VENT  OR  THROUGH  A 

C SIMPLE  FAN -RES I STANCE  VENTILATION  SYSTEM. 

C 

C***  SUBROUTINE  ARGUMENTS 
C 

C INPUT 

C 

C DIRS12 
C 

C YSLAB1 
C YSLAB2 
C XMSLAB 
C TSLAB 
C NSLAB 
C TU 
C TL 
C YFLOR 
C TEPS 
C 

C YCEIL 
C YLAY 
C QSLAB 
C RLAM 
C 

C PSLAB 
C MXPRD 
C NPROD 
C MXSLAB 
C 

C OUTPUT 

C 

C QSLVNT 
C 

C QSUVNT 
C 

C UFLW2 ( I , 1 , J ) , 1-1  OR  2,  J-l  OR  2 - MASS  FLOW  RATE  TO  UPPER 

(J-2)  OR  LOWER  (J-l)  LAYER  OF  ROOM  I DUE  TO  ALL  SLAB 
FLOWS  OF  VENT  [KG/S] 

UFLW2 (I , 2 , J ) , 1-1  OR  2,  J-l  OR  2 - ENTHALPY  FLOW  RATE  TO  UPPER 
C (J-2)  OR  LOWER  (J-l)  LAYER  OF  ROOM  I DUE  TO  ALL  SLAB 

C FLOWS  OF  VENT  [W] 


- A MEASURE  OF  THE  DIRECTION  OF  THE  ROOM  1 TO  ROOM 

2 FLOW  IN  EACH  SLAB 

- SLAB  HEIGHTS  IN  ROOM  1 ABOVE  DATUM  ELEVATION  [M] 

- SLAB  HEIGHTS  IN  ROOM  1 ABOVE  DATUM  ELEVATION  [M] 

- MASS  FLOW  RATE  IN  SLABS  [KG/S] 

- ABSOLUTE  TEMPERATURE  OF  THE  FLOW  IN  SLABS  [K] 

- NUMBER  OF  SLABS  BETWEEN  BOTTOM  AND  TOP  OF  VENT 

- UPPER  LAYER  TEMPERATURE  IN  EACH  ROOM  [K] 

- LOWER  LAYER  TEMPERATURE  IN  EACH  ROOM  [K] 

- HEIGHT  OF  FLOOR  IN  EACH  ROOM  ABOVE  DATUM  ELEVATION  [M] 

- MINIMUM  TEMPERATURE  DIFFERENCE  REQUIRED  TO  INITIATE 

GROWTH  OF  A NEW  LAYER  [K] 

- HEIGHT  OF  CEILING  IN  EACH  ROOM  ABOVE  DATUM  ELEVATION  [M] 

- HEIGHT  OF  LAYER  IN  EACH  ROOM  ABOVE  DATUM  ELEVATION  [M] 

- ENTHALPY  FLOW  RATE  IN  EACH  SLAB  [W] 

- FRACTION  OF  SOURCE  SLABS'  ENTHALPY  INFLOW  LOST 

TO  BOUNDING  WALLS 

- FLOW  RATE  OF  PRODUCT  IN  EACH  SLAB  [(UNIT  OF  PRODUCT/S] 

- MAXIMUM  NUMBER  OF  PRODUCTS  CURRENTLY  AVAILABLE. 

- NUMBER  OF  PRODUCTS 

- MAXIMUM  NUMBER  OF  SLABS  CURRENTLY  AVAILABLE. 


- TOTAL  RATE  OF  HEAT  TRANSFER  TO  LOWER  BOUNDING  SURFACES 
DUE  TO  ALL  SLAB  FLOWS  OF  VENT  [W] 

- TOTAL  RATE  OF  HEAT  TRANSFER  TO  UPPER  BOUNDING  SURFACES 
DUE  TO  ALL  SLAB  FLOWS  OF  VENT  [W] 
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C UFLW2 ( I , 3 , J ) , 1=1  OR  2,  J-l  OR  2 - OXYGEN  FLOW  RATE  TO  UPPER  (J-2) 

C OR  LOWER  (J=l)  LAYER  OF  ROOM  I DUE  TO  ALL  SLAB  FLOWS  OF 

C VENT  [(KG  OXYGEN) /S] 

C UFLW2 ( I , 3+K , J ) , 1=1  OR  2,  K=2  TO  NPROD , J=1  OR  2 - PRODUCT  K FLOW 

C RATE  TO  UPPER  (J=2)  OR  LOWER  (J=l)  LAYER  OF  ROOM  I DUE 

C TO  ALL  SLAB  FLOWS  OF  VENT  [(UNIT  PRODUCT  K)/S] 

C*END 

C 

C***  DEFINE  INDEXES  FOR  USE  WITH  PRODUCT  ARRAYS 
C 

INTEGER  L,U,M,Q,0, P 

PARAMETER  (L=l ,U=2 ,M=1 ,Q=2 ,0=3 , P=3) 

C 

INTEGER  DIRS12 (*) 

DIMENSION  YSLAB1(*) ,YSLAB2(*) ,XMSLAB(*) ,TSLAB(*) 

DIMENSION  QSLAB(*) 

DIMENSION  TU(*) ,TL(*) , YFLOR(*) ,YCEIL(*) , YLAY(*) 

DIMENSION  PSLAB(MXSLAB , *) 

DIMENSION  QSUVNT(*) ,QSLVNT(*) 

DIMENSION  UFLW2 ( 2 , MXPRD+2 , 2 ) , FF ( 2 ) 

C 

C***  INITIALIZE  OUTPUTS 
C 

DO  10  I = 1,  2 

DO  15  IPROD  = 1,  NPROD+2 
UFLW2 (I , I PROD , L)  = 0 . ODO 
UFLW2 (I , IPROD ,U)  = 0 . ODO 
15  CONTINUE 

QSUVNT ( I ) = O.ODO 
QSLVNT(I)  = O.ODO 
10  CONTINUE 
C 

C***  PUT  EACH  SLAB  FLOW  INTO  APPROPRIATE  LAYER  OF  ROOM  ITO 
C AND  TAKE  SLAB  FLOW  OUT  OF  APPROPRIATE  LAYER  OF  ROOM  I FROM 
C 

DO  20  N = 1,  NSLAB 
C 

C***  DETERMINE  WHAT  ROOM  FLOW  IS  COMING  FROM 
C 

IF (DIRS 12 (N) .EQ.1)THEN 
I FROM  = 1 
ITO  = 2 

ELSE  IF(DIRS12(N) ,EQ.O)THEN 
C 

C***  NO  FLOW  IN  THIS  SLAB  SO  WE  CAN  SKIP  IT 
C 

GO  TO  20 

ELSE  IF(DIRS12(N) . EQ. -1)THEN 
I FROM  = 2 
ITO  = 1 
ENDIF 
C 
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c***  APPORTION  FLOW  BETWEEN  LAYERS  OF  DESTINATION  ROOM 

C 

YLAYER  - YLAY(ITO) 

YLOW  - YFLOR(ITO) 

YUP  - YCEIL(ITO) 

TTS  - TSLAB(N) 

TTU  - TU(ITO) 

TTL  - TL(ITO) 

IF(ITO. EQ. 1)THEN 

YSLABT  - YSLAB1 (N) 

YSLABF  - YSLAB2 (N) 

ELSE 

YSLABT  = YSLAB2(N) 

YSLABF  = YSLABl(N) 

END  IF 
C 

C***  NO  UPPER  LAYER 
C 

IF (YLAYER . GE . YUP ) THEN 

IF(TTS . GT . TTL+1 . 00D0 ) THEN 
FU  - 1.0D0 
TTR  «=  TTL 
ELSE 

FU  - 0.0D0 
TTR  *=  TTL 
END  IF 
GO  TO  100 
ENDIF 
C 

C***  NO  LOWER  LAYER 
C 

IF (YLAYER . LE . YLOW) THEN 

IF (TTS . GE . TTU+1 . 00D0)THEN 
FU  - 1.0D0 
TTR  - TTU 
ELSE 

FU  - 0.0D0 
TTR  “ TTU 
ENDIF 
GO  TO  100 
ENDIF 
C 

C***  UPPER  LAYER  TEMP  > LOWER  TEMP 
C 

IF (TTU. GT. TTL) THEN 
C 

C***  DETERMINE  FU 
C 

IF(TTS.GE.TTU)THEN 
FU  - 1.0D0 

ELSE  IF (TTS. LE. TTL) THEN 
FU  - 0.0D0 
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ELSE 

FU  - (TTS  - TTL)/(TTU  - TTL) 

ENDIF 

C 

C***  DETERMINE  TTR 
C 

I F ( YSLABT . GT . YLAYER) THEN 
TTR  - TTU 

ELSE  IF(YSLABT . LT .YLAYER) THEN 
TTR  •=  TTL 
ELSE 

TTR  -=  (TTU+TTL) /2 . ODO 
ENDIF 
GO  TO  100 
C 

C***  UPPER  LAYER  TEMP  -=  LOWER  TEMP 
C 

ELSE  IF (TTU . LE . TTL) THEN 
IF(TTS . GT . TTU) THEN 
FU  = 1.0D0 
TTR  = TTU 

ELSE  IF(TTS . LT . TTU)THEN 
FU  = 0.0D0 
TTR  = TTL 
ELSE 

IF(YSLABT . GT .YLAYER) THEN 
FU  ■=  1.0D0 
TTR  - TTU 

ELSE  IF (YSLABT.LT. YLAYER) THEN 
FU  - 0.0D0 
TTR  - TTL 
ELSE 

FU  - . 50D0 
TTR  - TTL 
ENDIF 
GO  TO  100 
ENDIF 
ENDIF 

100  CONTINUE 

FL  - 1 . ODO-  FU 
FF(L)  - FL 
FF(U)  - FU 

FQ  - RLAM* ( 1 . ODO  - TTR/TTS ) 

FS  - 1.0D0  -FQ 
XMTERM  - XMSLAB(N) 

QTERM  - QSLAB(N) 

QSUVNT(ITO)  - QSUVNT(ITO)  + FU*QTERM*FQ 
QSLVNT(ITO)  - QSLVNT(ITO)  + FL*QTERM*FQ 
DO  110  I LAY  -1,2 

UFLW2(ITO,M,ILAY)  - UFLW2(ITO,M, ILAY)  + FF(ILAY)*XMTERM 
UFLW2 ( ITO , Q , ILAY)  - UFLW2 (ITO ,Q , ILAY)  + FF(ILAY)*QTERM*FS 
DO  120  IPROD  - 1,  NPROD 
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IP  - P + IPROD  - 1 

UFLW2(ITO, IP, ILAY)  - UFLW2 (ITO , IP , ILAY) 

1 + FF(ILAY)*PSLAB(N, IPROD) 

120  CONTINUE 
110  CONTINUE 
C 

C***  CASE  2A  FLOW  ENTIRELY  FROM  UPPER  LAYER  OF  ROOM  I FROM 
C 

IF ( YSLABF . GE . YLAY ( IFROM) ) THEN 

UFLW2 (IFROM,M,U)  - UFLW2 (IFROM, M,U)  - XMTERM 
UFLW2( IFROM, Q,U)  - UFLW2( IFROM, Q,U)  - QTERM 
DO  130  IPROD  ■=  1,  NPROD 
IP  «=  P + IPROD  - 1 

UFLW2 (IFROM, IP, U)  = UFLW2 (IFROM , IP ,U)  - PSLAB(N, IPROD) 
130  CONTINUE 

ELSE 
C 

C***  CASE  3A  FLOW  ENTIRELY  FROM  LOWER  LAYER  OF  ROOM  IFROM 
C 

UFLW2 (IFROM , M , L)  = UFLW2 ( IFROM ,M , L)  - XMTERM 
UFLW2(IFROM,Q,L)  = UFLW2(IFROM,Q,L)  - QTERM 
DO  140  IPROD  *=  1,  NPROD 
IP  = P + IPROD  - 1 

UFLW2  (IFROM,  IP  , L)  *=  UFLW2  (IFROM , IP  , L)  - PSLAB(N , IPROD) 
140  CONTINUE 

END  IF 
20  CONTINUE 
RETURN 
END 
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PLUGO 


Calculation  of  Mass,  Enthalpy,  Oxygen  and  Other  Product  of 
Combustion  Flows  to  Upper  and  Lower  Layers  Due  to  a Fire -Driven 
Plume 


DESCRIPTION 


Consider  an  instant  of  time  during  a two-layer  zone-type  simulation  of  a room 
fire  environment.  This  algorithm  calculates  the  rates  of  flow  of  mass, 
enthalpy,  oxygen  and  other  products  of  combustion  to  the  upper  and  lower 
layer  of  the  room  due  directly  to  a fire  in  the  room  and  its  fire -driven  plume 
[1],  The  fire  room  is  depicted  in  Figure  1. 

The  total  energy-release  rate  of  the  fire,  Q,  is  specified.  It  is  assumed 
that  the  combustion  zone  and  the  plume  radiate  at  an  effective  rate  ARQ,  where 
Ar  is  specified,  and  that 


^source  (1  ' 


(1) 


is  the  part  of  the  fire's  energy- release -rate  which  drives  the  convective 
buoyant  fire  plume  upward  toward  the  ceiling. 

The  fire  is  assumed  to  be  a source  of  oxygen,  02 , and  of  (NpR0D  - 1)  other 
user- specified  products  of  combustion  which  are  taken  account  of  in  the 
simulation.  Following  the  rule  of  [2]  it  is  assumed  that  the  rate  of  02 
generation  consumed  in  the  oxidation  process,  P02  source’  negative  ( i - e . , 
oxygen  is  consumed)  and  directly  proportional  to  Q. 

rate  of  02  generation  = P02> source  = 'F,Q  (2) 

F - 0 . 076(10" 6 ) [(kg  02)/(W-s)]  (3) 


Note  that  oxygenated  combustibles  would  have  additional  contributions  to  the 
rate  of  02  generation. 

It  is  assumed  that  the  other  products  of  combustion  are  generated  at  specified 
rates  PK  source > K = 2 to  NpR0D  [(unit  of  product  K)/s] . Typically,  the  rate 
at  which  mass  is  generated  in  the  combustion  process  is  negligible  compared  to 
other  mass  flow  rates  of  interest  in  the  simulation,  and  this  is  assumed  to  be 
the  case  here. 

If  the  lower  layer  is  of  non- zero  thickness  and  if  the  fire  is  in  the  lower 
layer,  then  it  is  assumed  that  a complete  description  of  the  lower- layer  room 
environment  at  the  instant  of  time  of  interest  is  known.  This  includes  the 
density,  absolute  temperature,  concentration-by-mass  of  02  [(kg  of  02)/(kg 
layer)],  and  concentrations -by-mass  of  the  other  NpR0D  products  of  combustion 
[(unit  of  product  k)/(kg  layer)].  These  are  designated  by  pL , TL , c02L,  and 
CK  L ’ K = 2 to  NpR0D , respectively.  The  02  and  other  product  concentrations 
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would  be  required  only  when  these  concentrations  were  actually  being 
predicted. 

Refering  to  Figure  1,  the  elevation  of  the  interface  separating  the  upper  and 
lower  layers  is  yLAYER . The  fire  is  located  at  elevation,  yFIRE , which  is 
assumed  to  be  between  or  at  the  elevation  of  the  room's  floor  and  ceiling, 
yFL00R  and  ycEiL>  respectively.  Consistent  with  the  idealized  point-source 
plume  model  of  [3],  it  is  appropriate  to  choose  yFIRE  as  the  lowest  elevation 
above  the  base  of  the  combustion  zone  where  relatively- free , unrestricted, 
lateral  entrainment  of  the  local  environment  into  the  combustion  zone  or  plume 
is  possible. 

It  is  assumed  that  the  total  fraction  of  Q which  is  lost  to  the  bounding 
surfaces  of  the  room  by  all  modes  of  heat  transfer  is  given  by  ATQ,  where  AT 
is  specified.  The  remaining  portion  of  Q is  Qjotal > i • e • > 


Qtotal  ~ (1  ” 


(4) 


As  explained  in  [1],  except  certain  special  configurations  this  is  assumed  to 
be  deposited  entirely  in  the  upper  layer.  When  the  fire  is  in  the  upper  layer 
this  is  all  added  to  the  upper  layer  directly  and  there  is  no  other  enthalpy 
flow  due  to  the  fire.  When  the  fire  is  in  the  lower  layer,  the  net  rate  of 
enthalpy  flow  into  the  upper  layer  due  to  the  fire  is  QT0TAL  plus  the  rate -of - 
enthalpy  flow  associated  with  the  gases  entrained  into  the  plume  from  the 
lower  layer  and  deposited  by  the  plume  into  the  upper  layer 

The  algorithm  assumes  no  gas- to-room- surface  losses  of  products  of  combustion. 


OUTPUT 


Mu  [MJ 

Rate  of  mass  flow  to  the  upper  [lower]  layer  of  the 
room  due  to  the  fire -driven  buoyant  plume.  [kg/s] 


^02,U  ^02,i] 

Rate  of  02  mass  flow  to  the  upper  [lower]  layer  of  the 
room  due  to  the  fire-driven  buoyant  plume.  [(kg  of 
02)/s] 

, U , L 5 > K = ^ t0  Mprod 

Rate  of  flow  of  product  of  combustion  K to  the  upper 
[lower]  layer  of  the  room  due  to  the  fire -driven 
buoyant  plume.  [(unit  of  product  K)/s] 
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Qu  [QJ 

Rate  of  enthalpy  flow,  based  on  absolute  temperature, 
to  the  upper  [lower]  layer  of  the  room  due  to  the 
fire -driven  buoyant  plume.  [W] 


Qtotal 

Net  rate  of  enthalpy  flow  to  the  combined  upper  and 
lower  layers,  (1  - AT)Q.  [W] 


INPUT 


Specific  heat  at  constant  pressure  of  the  layers 
(suggest  using  103  W-s/(kg-K)  which  is  Cp  of  air). 
[W-s/(kg-K)] 


C0  2 , L 


Concentration  of  02  in  the  lower  layer.  [(kg  of 
02)/(kg  of  layer)] 


2 to  Nr 


Concentration  of  product  of  combustion  K in  the  lower 
layer.  [(unit  of  product  K)/(kg  of  layer)] 


Number  of  products  of  combustion,  including  oxygen, 
being  tracked  in  the  simulation. 


^K, SOURCE’  K »=  2 to  NpR0D 

Rate  of  product  of  combustion  K generated  by  the  fire 
source.  [(unit  of  product  K)/s] 


Q 


Total  rate  of  energy  release  of  the  fire.  [W] 
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T, 


Absolute  temperature  of  the  lower  layer.  [K] 


yCEILING/yLAYER/yFIRE 

Elevation  of  the  ceiling/layer/fire  above  the  datum 
elevation.  [m] 


A small  temperature  difference,  between  the  average 
plume  temperature  at  the  ceiling  and  TL , below  which 
it  is  assumed  that  an  upper  layer  will  not  grow  from 
zero  thickness.  [K] 


Fraction  of  Q radiated  by  the  combustion  zone  and 
plume,  0 < Ar  <1. 

At 

Total  fraction  of  Q lost  to  the  bounding  surfaces  of 
the  room  by  all  modes  of  heat  transfer,  0 < AT  <1. 


density  of  the  lower  layer.  [kg/m3] 


CALCULATION 


The  calculations  are  presented  in  the  flow  diagram  of  Figure  2 . 
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SUBROUTINE  VARIABLES 


Nomenclature  used  above 
subroutine : 

correspond 

to  the  following  nomenclature  used 

Cp 

- 

CP  [W- s/ (kg-K) ] 

FACTOR 

- 

0 . 076(10" 6 ) [(kg  of  02)/(W-s)] 

g 

- 

9 . 8 [m/s2 ] 

Ml  [MJ 

- 

ML,  MU  [kg/s] 

M 

1 iP  L U M E 

- 

MPLUM  [kg/s] 

, SOURCE 

- 

PSOURC(K)  [(unit  of  product  K)/s] 

Ql>  Qu 

- 

QL,  QU  [W] 

^SOURCE 

- 

QSOURC  [W] 

Qtoxal 

- 

QTOT  [W] . 

- 

TL  [K] 

yCEILING 

- 

YCEIL  [m] 

YV  I RE 

- 

YFIRE  [m] 

Flayer 

- 

YLAY  [m] 

Z 

- 

Z [m] 

^^E  P S 

- 

DTEPS  [K] 

atplume 

( 

- 

DT  [K] 

Pl 

- 

RHOL  [kg/m3 ] 

Ar,  at 

- 

ALAMR,  ALAMT 

PREPARED  BY 

Leonard  Y.  Cooper  and  Glenn  P.  Forney 
March  1988 


PLUGO  - 5 


/V7\  vssssssss/ss 


Figure  1.  Fire  environment  in  a room  of  fire  origin 
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Figure  2.  Flowchart  of  PLUGO  calculations 
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SUBROUTINE  PLUGO 


SUBROUTINE  PLUGO ( 

I CP , CL , NPROD , PSOURC , Q , TL , TU , 

I YCEIL , YLAY , YFIRE , DTEPS , ALAMR , ALAMT , RHOL , 

0 XMU , XML , PU , PL , QU , QL , QTOT ) 

IMPLICIT  DOUBLE  PRECISION  (A-H.O-Z) 

C*BEG 

C 

C***  PLUGO  GENERIC  - ROUTINE  TO  CALCULATE  THE  MASS,  ENTHALPY  AND  PRODUCT 
C FLOWS  DUE  TO  ONE  PARTICULAR  PLUME. 

C 

C***  SUBROUTINE  ARGUMENTS 
C 

C INPUT 
C 

C CP  SPECIFIC  HEAT  [W*S/(KG*K) ] 

C CL  CONCENTRATION  OF  EACH  PRODUCT  IN  LOWER  LAYER 
C [UNIT  OF  PRODUCT/ (KG  LAYER)] 

C NPROD  NUMBER  OF  PRODUCTS  OF  COMBUSTION 

C PSOURC  AMOUNT  OF  PRODUCT  RELEASED  (OR  ABSORBED)  BY  THE  FIRE 
C Q ENERGY  RELEASE  RATE  OF  FIRE  [W] 

C TL  TEMPERATURE  OF  LOWER  LAYER 

C TU  TEMPERATURE  OF  UPPER  LAYER 

C YCEIL  HEIGHT  OF  CEILING  ABOVE  DATUM  ELEVATION  [M] 

C YLAY  HEIGHT  OF  LAYER  ABOVE  DATUM  ELEVATION  [M] 

C YFIRE  HEIGHT  OF  FIRE  ABOVE  DATUM  ELEVATION  [M] 

C DTEPS  MINIMUM  TEMPERATURE  OF  FLOW  THAT  WILL  CAUSE  FORMATION 
C OF  AN  UPPER  LAYER  [K] 

C ALAMR  FRACTION  OF  Q RADIATED  BY  THE  COMBUSTION  ZONE  AND  PLUME 

C ALAMT  TOTAL  FRACTION  OF  Q LOST  TO  THE  BOUNDING  SURFACES  OF  THE  ROOM 

C BY  ALL  MODES  OF  HEAT  TRANSFER 

C RHOL  DENSITY  OF  LOWER  LAYER  [KG/M**3] 

C 

C OUTPUT 
C 

C XMU  MASS  FLOW  RATE  INTO  UPPER  LAYER 

C XML  MASS  FLOW  RATE  INTO  LOWER  LAYER 

C PU  PRODUCT  FLOW  RATE  INTO  UPPER  LAYER 

C PL  PRODUCT  FLOW  RATE  INTO  LOWER  LAYER 

C QU  ENTHALPY  FLOW  RATE  INTO  UPPER  LAYER 

C QL  ENTHALPY  FLOW  RATE  INTO  LOWER  LAYER 

C QTOT  TOTAL  ENTHALPY  FLOW  RATE 

C 

C*END 

PARAMETER  (G-9 . 8D0 , FACTOR= . 076D-6) 

DIMENSION  PU(*),  PL(*),  PSOURC (*) , CL(*) 

C 

QTOT  - (1.0D0  - ALAMT ) *Q 
QSRC  = (1.0D0  - ALAMR) *Q 
Z - YLAY  - YFIRE 

IF(NPROD . GE . 1 ) PSOURC ( 1)  - -FACTOR*Q 
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IF( Z . LT . 0 . ODO . OR . ( Z . EQ . 0 . ODO . AND . YFIRE . LT . YCEIL) ) THEN 
C 

C FIRE  IN  UPPER  LAYER;  NO  PLUME  CALCULATIONS  REQUIRED; 

C NO  CHANGE  TO  LOWER  LAYER 
C 

XML  = 0.0D0 
QL  = 0.0D0 
XMU  = O.ODO 
DO  10  IPROD  = 1,  NPROD 
PL ( I PROD)  = O.ODO 
PU (IPROD)  = PSOURC( IPROD) 

10  CONTINUE 
QU  = QTOT 
RETURN 
END  IF 

IF ( Z.EQ. O.ODO. AND. YFIRE. GE. YCEIL) GO  TO  300 
C 

C FIRE  IN  LOWER  LAYER;  AND  BELOW  INTERFACE; 

C USE  PLUME  CALCULATIONS  FOR  XMPLUM,  TPLUME 
C 

GZ  = (G*Z**5)** . 5D0 

QSTAR  = QSRC/(RHOL*TL*CP*GZ) 

XMPLUM  = .21D0*RHOL*GZ*QSTAR**(1.0D0/3.0D0) 

IF (YCEIL . LE . YLAY . AND . XMPLUM . NE . 0 . ODO ) THEN 
DTPLUM  = QSRC/ ( CP*XMPLUM) 

IF(DTPLUM. LT . DTEPS)GO  TO  300 
END  IF 

XMU  - XMPLUM 
DO  20  K - 1,  NPROD 

PU(K)  = CL ( K ) *XMPLUM  + PSOURC(K) 

PL(K)  = - CL ( K ) *XMPLUM 
20  CONTINUE 

QU  - CP*XMPLUM*TL  + QTOT 
XML  - -XMPLUM 
QL  = - CP*XMPLUM*TL 
RETURN 
300  CONTINUE 
C 

C NO  CHANGE  TO  UPPER  LAYER  SINCE  FIRE  IS  IN  LOWER  LAYER  AND: 

C 1)  IT  IS  IMMEDIATELY  BELOW  BOTH  CEILNG  AND  INTERFACE,  OR 

C 2)  INTERFACE  IS  AT  CEILING  AND  PLUME  IS  TOO  WEAK  TO  START  NEW  LAYER 

C 

XMU  - O.ODO 
QU  = O.ODO 
DO  30  K = 1,  NPROD 
PL(K)  = PSOURC(K) 

PU(K)  = O.ODO 
30  CONTINUE 
XML  - O.ODO 
QL  = QTOT 
RETURN 
END 
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VENTF 


Calculate  Flow  Through  a Fan/Duct  Forced  Ventilation  System 
Between  Two  Rooms 


DESCRIPTION 

Consider  an  instant  of  time  during  the  two -layer,  zone -type  simulation  of  a 
multi -room  compartment  fire  environment.  This  algorithm,  VENTF,  calculates 
the  direction  and  rate  of  flow  of  mass,  enthalpy,  oxygen  and  products  of 
combustion  through  a fan/duct  forced  ventilation  system  connecting  any  two 
rooms,  or  a room  and  the  outside  environment  [1].  The  calculation  uses  a 
separate  algorithm,  FANRES  [2],  which  performs  the  system  flow  simulation  in 
cases  where  the  flow  rate  is  not  specified. 

Designate  the  two  rooms  or  the  room  and  the  outside  environment  as  rooms  1 and 
2 and  refer  to  the  Figure  1.  Define  room  1 [room  2]  as  the  room  connected  to 
the  system  inlet  [outlet].  Here,  the  system  inlet  [outlet]  is  defined  as  the 
side  of  the  fan/duct  flow  system  which  corresponds  to  the  flow  inlet  [outlet] 
under  normal  fan  operating  conditions.  Normal  operating  conditions  lead  to 
what  will  be  referred  to  here  as  flow  in  the  normal  or  forward  direction. 

The  elevation  of  the  inlet  and  outlet  end-point  of  the  fan/duct  system  are 

characterized  by  the  elevation  of  their  midpoints,  ySYSTEM  i>  ^-n  roonl  I = 1 or 

2,  above  a datum  elevation.  The  vent  openings  at  the  duct  end-points  are 
assumed  to  be  rectangular  with  sides  aligned  vertically  and  horizontally.  Hj , 
the  depth  of  these  vents  in  room  i,  are  also  specified.  These  are  required  in 
order  to  determine  the  effect  on  the  fire  environment  of  the  system  flow  at 
instants  of  time  when  flow  through  the  system  is  drawn  partially  from,  or 
supplied  partially  to  both  the  upper  and  lower  layer  of  room  I. 

A complete  characterization  of  the  instantaneous  environment  in  each  of  the 
two  rooms  is  provided  as  input.  Using  this,  the  algorithm  first  computes  the 
densities  of  the  flow  in  the  ventilation  system  duct  near  its  entry  region 
under  the  conditions  of  either  forward  flow  (flow  enters  at  the  designed 

system  inlet)  or  backward  flow  (flow  enters  at  the  designed  system  outlet) . 

Using  these  two  densities,  certain  specified  characteristics  of  the  fan/ duct 
system,  and  FANRES  [2] , the  algorithm  estimates  the  instantaneous  direction  of 
flow  through  the  system,  indicated  by  the  parameter  IDIR,  and  the  mass-flow 
rate,  MDUCT . Assuming  an  adiabatic  duct  flow,  the  algorithm  then  calculates: 
the  temperature  of  the  duct  flow,  Tduct , and  the  rates  of  flow  through  the 
duct  of  enthalpy,  oxygen,  and  NpR0D  - 1 other  products  of  combustion,  QDUCT , 
Pduct,02>  *duct,k>  K - 2 to  NpR0D , respectively.  Finally,  the  algorithm 
computes:  the  number  and  rates  of  flow  of  uniform-property  flow  slabs  leaving 
the  source  room  (one  or  two  flow  slabs  depending  on  whether  the  flow  is 
extracted  from  one  or  both  of  the  layers)  and  the  number  and  rates  of  flow  of 
identically-uniform-property  flow  slabs  entering  the  recieving  room.  If  there 
is  only  one  of  these  latter  slabs  entering  the  recieving  room,  then  it  enters 
the  upper  or  lower  layer.  If  there  are  two  such  flow  slabs,  then  one  enters 
the  upper  layer  and  the  other  the  lower  layer. 

The  specified  characteristics  of  the  fan/duct  system  include  pairs  of  values, 
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[VF  Jt  ApF  j ] , J - 1 to  JTERM,  which  identify  points  on  the  fan  operating 
curve . 

If  JTERM  - 1,  then  the  single  VF  z input  is  a specified  volume -flow- rate . 

Under  these  circumstances,  the  corresponding  ApF  z input  is  not  required  and 
FANRES  in  not  used. 

If  JTERM  > 1,  then  the  specified  [VF  Jf  ApF  j]  values  indicate  the  volume- 
flow-rate  delivery  of  the  fan,  Vp  j,  at  discrete  cross-fan  pressure 
differences,  ApF  j.  Referring  to  Figure  2,  ApF  j is  defined  as  the  cross-fan 
pressure  difference,  ApF , at  point  J on  the  fan’ curve  where 

aPf  ” 5Pf,2  * 5Pf, i C1) 

and  where  5pF  z and  5pF  2 are  the  pressures  that  would  be  measured  at  the  fan 
inlet  and  outlet  sections,  respectively,  relative  to  the  datum  pressure.  In 
specifying  the  [VF  j,  ApF  j]  data,  the  ApF  j are  required  to  be  monotonically 
decreasing  with  increasing  VF  j . A sketch’ of  such  a generic  fan  operating 
curve  is  presented  in  Figure  2.  In  the  calculations,  the  fan  curve  is 
approximated  by  linear  interpolation  between  and  linear  extrapolation  beyond 
the  fan-curve  data  points. 

If  JTERM  > 1,  then  effective  flow  resistances  of  the  duct  system,  Rx  and  R2 
for  forward  and  backward  flow,  respectively,  must  be  specified.  (The 
definitions  for  the  Rj ' s are  presented  in  the  list  of  inputs.)  Finally,  the 
midpoint  elevation  of  the  two  vents  at  the  system  endpoints  are  specified. 
Further  discussion  on  the  fan/ duct  system  characteristics  and  the  basis  of  the 
FANRES  algorithm  [2]  are  presented  in  [3]. 

As  mentioned  above,  flow  through  the  system  is  determined  by  the  instantaneous 
environments  local  to  the  vents  at  the  two  system  endpoints.  As  depicted  in 
Figure  1 it  is  assumed  that  these  two  environments  are  described  within  the 
context  of  a two- layer  model  of  the  rooms'  fire  environments.  In  each  of  the 
rooms  I,  I - 1 or  2,  VENTF  computes  the  local  environment  from  specified 
values  of:  pressure  at  the  floor  elevation,  $pFL00R  j . relative  to  a datum 
pressure;  elevations  of  the  floor  and  smoke  layer  interface,  yFLOoR  i ancl 
yLAYER  i * above  a datum  elevation;  absolute  temperature  and  density  of  the 
upper  and  lower  layer,  Ty  z , TL  z and  z , pL  z , respectively;  and 
concentrations  of  oxygen  and  products  of  combustion  in  the  upper  and  lower 
layer,  cU02I,  cL  02>I  and  c^ ( R ( j , cL  K z , K - 2 to  NpR0D , respectively. 

As  discussed  in  [1]  and  [2],  at  some  instants  of  time  in  a fire  scenario  the 

existance  of  two  solutions  to  the  fan/duct  flow  equations,  a forward-  and 
backward-flow  solution,  is  possible.  However,  at  a particular  instant  of  time 
in  a simulation,  it  not  known  apriori  whether  or  not  there  is  a single  or 
multiple  solution.  For  this  reason  a final  input  variable , ID IR0LD , is 
required  which  identifies,  in  the  event  of  a multiple  solution,  the  preferred 
direction  of  the  flow.  When  carrying  out  simulations  of  time -dependent  fire 
environments  it  is  reasonable  to  choose  the  preferred  flow  direction  for  each 
fan/duct  system  as  the  flow  direction  computed  at  the  previously  simulated 

instant  of  time  (or,  as  the  normal  flow  direction,  if  the  previous  flow  rate 

is  zero  identically) . 
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Characterizing  the  Entering  Duct  Flow  When  the  Elevation  of  the  Layer 

Interface  is  Between  the  Top  and  Bottom  of  the  Vent  Opening 


When  the  layer  interface  is  not  between  the  top  and  bottom  of  the  vent  opening 
the  properties  of  the  entering  duct  flow  are  identical  to  the  specified 
properties  of  the  lower  layer  (if  yLAYERiIFR  ^ Ysystem.ifr  + hifr/2>  or  uPPer 
layer  (if  yLAYER>IFR  < Ysystem.ifr  " hifr/2>  where  IFR-  standing  for  I FROM,  is 
the  number  (1  or  2)  of  the  source  room  from  which  the  flow  enters  the 
ventilation  system.  (ITO  will  be  used  below  as  the  number  of  the  recieving 
room  into  which  the  ventilation  flow  is  deposited.)  When  the  layer  interface 
is  between  the  top  and  bottom  of  the  vent  opening  special  considerations  are 
required. 

Figure  3 is  an  illustration  of  the  flow  entering  the  duct  in  room  IFR  when  the 
layer  interface  is  between  the  top  and  bottom  of  the  vent  opening.  Define 
5pDUCT  as  the  pressure  relative  to  the  datum  pressure  at  a location  somewhat 
inside  the  duct,  upstream  of  significant  duct- flow  development  and  mixing 
between  the  upper-  and  lower-layer  flow  contributions. 


Assume  that  the  flow  through  the  vent  and  into  the  entry  region  of  the  duct  is 
similar  to  the  flow  through  an  orifice  or  nozzle.  Assume  also  that  HIFR  is 
small  enough  so  that  the  absolute  value  of  the  hydrostatic  pressure  difference 
over  the  height  of  the  duct  is  small  compared  to  the  absolute  value  of  the 
pressure  drop,  |ApENTRY|  ” |5pDUCT  ' ^Psystem  ifrI»  across  the  duct  entrance 
region.  Then  the  velocity  of  the  upper  and  lower  layer  flow- slab  components 
of  the  flow  in  the  entrance  region  of  the  duct  are  both  proportional  to 
| ApENTRY 1 1 / 2 • In  particular,  the  mass  - flow-rate  of  the  two  flow  slabs  from 
the  upper  and  lower  layer,  MgLABU  IFR  and  MsLABL  iFR»  respectively,  are  given 

by 

^slabu.ifr  “ bCDbu , IFR (I^Pentry IPu , ifr)1/2 

(2) 

^SLABL.IFR  “ bCDhL > IFR ( |ApENTRY |pL > IFR)1/2 

In  the  above,  b is  the  width  of  the  vent  and  CD  is  the  flow  coefficient. 

Also,  as  illustrated  in  Figure  3 , hu  x [hL  j]  is  the  absolute  value  of  the 
difference  in  elevation  in  room  I between  the  top  [bottom]  of  the  vent  and  the 
layer  interface. 


bu , i “ Ysystem.i  + Hi/2  ‘ Ylayer.i 

(3) 

^L.I  HI  bU(I  - yLAYER,I  " (Ysystem.i  ■ Hl/2) 

From  Eq.  (2)  it  follows  that 

^SLABU  , IFr/^SLABL  . IFR  ” (bU  , I FR  /hL  ,1  F R ) ( P\j  , i F R / P\.  , j F R ) 1 ' 2 (^) 

Assume  that  once  in  the  duct  the  two  streams  mix  adiabatically  to  form  a 
uniform-property  duct  flow  of  temperature,  Tduct , and  density,  PDUCT . where 
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the  upper-  and  lower- layer  inlet  flows  can  all  be  modeled  as  the  same  perfect 
gas , i . e . , 

/5DUCT^,DUCT”^U,IFR,^U(IFR‘=^L,IFR^'l,IFR  ( 5 ) 

Adiabatic  mixing  requires 

^SLABU , IFRTU , IFR  + ^ L A B L , I F R TL , I F R “ IMDUCtITDUCT 

“ (^SLABU.IFR  + ^LABL  , IFR  )TDUCT 

(6) 


Now  Eqs . (4) -(6)  lead  to 

^DUCT  “ ^DUCT.IFR  = ^L.IFR^  + , I Fr/^L  , I FR  , I Fr/^L  , I FR  ) 1 ' 2 ^ / 

[1  + (hu ( IFR /hL , IFR ) , IFr/^U , IFR ^ ' 2 ] 

(7) 

where  pDUCT  IFR  is  the  density  of  the  duct  flow  when  the  source  room  is  room 
IFR.  Thus,’ if  the  flow  direction  (i.e.,  the  value  of  IFR)  is  known,  the  above 
model  for  flow  entry  into  the  duct  yields  an  Eq.-(7)  estimate  of  PDUCt»ifr 
without  specific  knowledge  of  the  value  of  MDUCT  or  its  entry  flow-slab 
components,  MsLABU  IFR  and  MsLABL  IFR.  This  result  will  be  used  below  in  the 
FANRES  algorithm  wfiich  requires  values  of  pD  IFR  under  possible  conditions  of 
both  forward  or  backward  flow  through  the  ventilation  system,  i.e.,  for  both 
possibilities  IFR  - 1 or  2 . 

With  the  solution  PDUCT  IFR  for  both  possibilities  IFR  - 1 and  2,  it  is 
possible  to  obtain  the  correct  solution  for  IFR  and  the  correct  associated 
value  for  the  flow  rate.  If  JTERM  > 1,  this  will  require  use  of  the 
algorithm/subroutine  FANRES  [2]. 

With  the  correct  value  of  IFR,  the  solution  to  pDUCT  is  determined  in  Eq.  (7). 
Then,  Eq.  (5)  leads  to  Tduct 

tduct  " P\i , ifrtu  , ifr/^duct 

Once  Mduct  is  available,  Eqs.  (4),  (6),  and  (8)  lead  to  MsLABU,ifr  and 

^SLABL , IFR 

^SLABL.IFR  “ I^DUCTl^DUCT/^U.IFR^ku.IFR  /^L  ,IFr)(^U,IFr/^L.IFr)1/2] 


(9) 
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^SLABU , IFR 


INDUCT  I " ^SLABL.IFR 

INDUCT  lTDUCT/^TL,IFRt  (hL  , IFr/^U  .IFR^^L.IFr/^U.IFR^^^ 


(10) 


The  Flow  in  the  Duct  and  the  Flow  Slabs  Entering  the  Recieving  Room 

The  above  ideas  lead  to  the  solution  to  the  duct -entry  flow  problem.  With 
this,  the  assumed  fully-mixed  duct-flow  condition  leads  first  to  a solution 
for  all  of  the  duct-flow  properties,  and  then  to  a solution  for  the  number  (1 
or  2)  and  properties  of  the  flow  slabs  which  enter  the  recieving  room. 


OUTPUT 

*DIR 

A measure  of  the  direction  of  the  room  1- to -room  2 
system  flow.  For  flow  from  room  1 to  room  2 (normal 
inlet  to  normal  outlet)  or  from  room  2 to  room  1 
(normal  outlet  to  normal  inlet,  i.e.,  a condition  of 
fan  backflow),  the  value  of  ID1R  is  1 or  -1, 
respectively.  IDIR  is  0 if  the  room-to-room  flow  is 
identically  zero. 

^SLABF,N>  N - 1 to  Nslabf  < 2 [MsLABT.N^  N - 1 to  Nslabt  < 2] 

Total  mass  flow  rate  in  slab  N leaving  the  source 
room,  room  IFR  [entering  the  recieving  room,  room 
ITO].  [kg/s] 

^duct 

Mass  flow  rate  through  the  duct.  [kg/s] 

nprod 


Number  of  products  of  combustion  being  tracked  in  the 
simulation. 
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Number  of  slabs  of  flow,  1 or  2,  leaving  the  source 
room,  room  IFR,  and  entering  the  system  duct.  If 
NSLabf  ” 1 then  the  flow  is  totally  from  the  upper  or 
lower  layer.  If  NSLABF  - 2 then  each  layer  supplies 
one  of  the  slabs. 


nslabt 

Number  of  slabs  of  flow,  1 or  2,  supplied  to  the 
recieving  room,  room  ITO,  and  leaving  the  system  duct 
If  Nslabt  *=  1 then  all  of  the  flow  enters  either  the 
upper  or  lower  layer.  If  NSLABT  — 2 then  one  of  the 
slabs  will  enter  the  upper  layer  and  the  other  slab 
will  enter  the  lower  layer. 

^duct,k>  K - 2 to  NpR0D 

Product  K flow  rate  through  the  duct.  [(unit  of 
product  K)/s] 


^duct , 02 


Mass  flow  rate  of  oxygen  flowing  through  the  duct, 
[(kg  of  oxygen) /s] 


LABF.K.N’  w uo  l’S  L A B F ~ ^ c LO  iNPR0D 

L ab T , k , N ’ N - 1 to  Nslabt  <2,  K - 2 to  NpR0D ] 

Product  K flow  rate  in  slab  N leaving  the  source  room 
room  IFR  [entering  the  recieving  room,  room  ITO] . 
[(unit  of  product  K)/s] 


**SLABF  . 02  , N ’ N - 1 to  NgLABp  < 2 [ PS  L AB  T , 0 2 , N > N - 1 to  NgLABT  < 2] 

Oxygen  flow  rate  in  slab  N leaving  the  source  room, 
room  IFR  [entering  the  recieving  room,  room  ITO] . 

[ (kg  of  02)/s] 

Qduct 


Enthalpy  flow  rate  through  the  duct  based  on  absolute 
temperature  Tduct . [W] 
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Qslabf.N’  N - 1 to  Nslabf  < 2 [QSlabt,n>  N - 1 to  NSLABT  < 2] 

Enthalpy  flow  rate  in  slab  N leaving  the  source  room, 
room  IFR,  based  on  the  absolute  temperature,  TSLABF  N 
[entering  the  recieving  room,  room  ITO,  based  on  the 
absolute  temperature,  Tslabt  „].  [W] 


tduct 


Absolute  temperature  of  the  flow  in  the  duct.  [K] 


^slabf,n>  N = 1 to  Nslabf  < 2 [TSLABF  N ; N - 1 to  NSLABF  < 2] 

Absolute  temperature  of  the  flow  in  slab  N leaving  the 
source  room,  room  IFR  [entering  the  recieving  room, 
room  ITO] . [K] 

L AB  F , N ’ N ” 1 t°  Nslabf  < 2 [ySLABT , N ’ N - 1 to  NSLABT  < 2] 

Elevations  above  the  datum  elevation  of  the  midpoints 
of  slab  N leaving  the  source  room,  room  IFR  [entering 
the  recieving  room,  room  ITO] . [m] 

INPUT 

ACC 

Absolute  error  tolerance  for  the  fan/duct  system's 
flow  rate  calculation.  This  parameter  is  not  used 
when  the  flow  rate  is  specified,  i.e.,  when  JTERM  = 1. 

cl  , k , i t cu , k , 1 1 ’ K - 2 to  NPR0D ; I - 1 or  2 

Concentration  of  product  of  combustion  K in  lower 
[upper]  layer  of  room  I.  [(unit  of  product  K)/(kg  of 
layer) ] 

CL  ,02,1  t CU , 0 2 , I ] ’ 1 “ ^ °r  ^ 

Concentration  of  oxygen  in  lower  [upper]  layer  of  room 
I.  [(kg  of  oxygen)/ (kg  of  layer)] 
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Specific  heat  at  constant  pressure  of  the  vent  flow. 
[W-s/(kg-K)]  (suggest  103  W-s/(kg-K)  for  air  as 
default) 


^•DIROLD 

A measure  of  the  preferred  direction  for  the  solution 
flow  in  the  event  of  multiple  flow  solutions,  i.e., 
the  preferred  value  for  IDIR. 

Hj ; I - 1 or  2 

Height  of  the  vent  in  room  I.  Note  that  Hj  must 
satisfy:  Hx  < 2(yCEILING  - ys y stem, i ) and  Hi  * 
2(ysYSTEM,I  ' yFLOOR , I ^ ‘ 

JTERM  > 1 

Number  of  pairs  of  input  data  [VF  j , ApF  3 ] used  to 
approximate  the  fan  curve . 

Nprod 

Number  of  products  of  combustion  being  tracked  in  the 
simulation. 


Ri  [Rj] 


(Required  if  JTERM  >1.)  Effective  resistance  of  the 
duct  for  forward  [backward]  flow.  If  MDUCT  z , is  the 
mass  flow  rate  through  the  fan/duct  system,  then  Rx  is 
defined  as 

RI  ” IAPf  + (^PsYSTEM.l  ' ^PsYSTEM , 2 ^ 

+ ^DUCT , I SCYsYSTEM , 1 " ySYSTEM,2^ |1/2/|MDUCT>I I 

where  App  is  defined  in  Eq.  (1),  and  ^PSystem,i  and 
5pSYSTEM  2 are  t^ie  pressures  in  rooms  1 and  2, 
relative  to  the  datum  pressure,  at  the  elevations, 

ySYSTEM.l  and  ySYSTEM , 2 [Pa1/2kg  XS  - 

{ kg/ (m • s^ ) } 1 / 2 kg” * s ] 
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Vi  [Tl>i1;  I - 1 or  2 


Absolute  temperature  of  the  upper  [lower]  layer  in 
room  I.  [K] 

Vj,  J - 

1 to  JTERM 

Volume- flow- rate  delivered  by  the  fan  when  cross -fan 
pressure  difference  is  at  the  specified  value  ApF  j . 

A VF  j value  is  positive  when  flow  is  in  the  normal 
direction  of  fan  operation.  The  pairs  of  (VF  j, 

ApF  j ) data  would  be  obtained  from  fan  operating 
characteristics  provided  typically  by  the  fan 
manufacturer.  The  ApF  j are  required  to  be 
monotonically  decreasing  with  increasing  VF  j . 

The  ApF  corresponding  to  a given  volume- flow- rate 
through  the  fan,  VF , is  calculated  by  the  FANRES 
algorithm/subroutine  [2].  This  is  done  by 
interpolating  between  the  (Vp  Jf  ApF  j)  data  points  or 
by  extrapolating  beyond  these  for  VF  < VF  1 and  VF  < 

Vp  j x erm  • If  JTERM  - 1,  then  the  Vp  x input  is  the 
specified  volume- flow- rate , the  corresponding  ApF  , 
input  is  not  required,  and  FANRES  is  not  used.  [nr/s] 

^FLOOR , I » 

I - 1 or  2 

Elevation  of  the  floor  in  room  I above  the  datum 
elevation.  [m] 

Flayer , i • 

I - 1 or  2 

Elevation  of  layer  in  room  I above  the  datum 
elevation.  [m] 

^SYSTEM , I 

; I - 1 or  2 

Elevation  of  the  middle  of  the  system  endpoint  in  room 
I above  the  datum  elevation.  [m] 

aPf.j>  j 1 

- 1 to  JTERM 

Values  of  the  cross -fan  pressure  difference  at  the 
JTERM  points  used  to  approximate  the  fan  operating 
curve  (volume  flow- rate  vs  cross -fan  pressure 
difference).  See  the  above  VF  j entry.  [Pa  - 
kg/(m-s2) ] 
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^Pfloor , i * I - 1 or  2 

Pressure,  above  the  datura  pressure,  at  the  floor  in 
room  I.  [Pa  = kg/(m-s2)] 


Pl.i  [Pu.i] ; I = 1 or  2 

Density  of  lower  [upper]  layer  in  room  I.  [kg/m3] 


CALCULATIONS 

1.  If  JTERM  - 1 and  VF  x =0,  then  there  is  no  forced  vent  flow  and  no 
adjustments  are  made  to  the  environments  in  rooms  1 and  2.  If  JTERM  = 1 
and  VF  x is  non- zero,  or  if  JTERM  > 1,  then  continue  with  calculation 
step  2. 

2.  For  I - 1 and  2: 

a.  Calculate  5psysTEM>I: 

i.  if  ySysTEM,i  > Flayer , i ■ 

5Psys  T EM , I “ PfLOOR.X  ' Pl , I S(yLAYER , I " yFLOOR.l) 

” P\s , i S(ys ystem  , i ~ Flayer  , i ) (10) 

ii'.  if  ySYSTEM,I  ^ yL AYER  , I * 

5Psys  TEM , I “ PFLOOR , I 

■ Pl . i S(Ysystem , I " ^FLOOR , I ) (^^) 

b.  Compute  h0  x and  hL  x from  Eq.  (2). 

c.  Although  the  flow  direction  is  not  yet  known,  the  density  in  the 
duct  near  the  room  IFR  flow  entry  is  computed  for  the  possible 
flow  condition  IFR  - I,  for  I - 1 or  2. 

1*  If  YSTEM,  I + Hl/2  ^ yLAYER,i: 

the  duct  flow  is  entirely  from  the  lower  layer  and 
Pduct , i “ Pl . i 

ii . if  ys y stem , i - Hi/2  ^ yLAY  ER  , I • 

the  duct  flow  is  entirely  from  the  upper  layer  and 
Pduct  , i “ P u . i 
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iii.  else: 


the  duct  flow  is  from  both  the  upper  and  lower  layers 
and  pDUCT  j is  computed  from  Eq.  (7) 

3.  Determine  MDUCT , IDIR,  IFR,  and  ITO.  If  JTERM  > 1,  then  use  the  FANRES 
algorithm.  If  JTERM  - 1,  then  MDUCT  is  obtained  from  the  specified 
volumetric  flow  rate,  VF  1# 

a.  if  JTERM  = 1: 

i.  if  VF  1 > 0: 

^DUCT  “ ^DUCT,  l^F,  1 > ^DIR  “ ^ ’ an<^  IFR  “ I 

ii.  if  VF>1  < 0: 

^DUCT  “ '^DUCI,2^F,1  ^DIR  = "i»  an<^  IFR  = ^ 

b.  if  JTERM  > 1 calculate: 

SYSTEM  ” ySYSTEM.l  * ySYSTEM,2 
^PSYSTEM  “ ^PSYSTEM.1  ’ ^PsYSTEM.2 
and  then: 

CALL  FANRES(M30LN  N,  Nsoln  , ACC,  JTERM,  R, , Vp  . J • ApF  a , 

^SYSTEM*  ^SYSTEM’  ^DUCT.I^ 

i-  if  Nsoln  «=  1,  then  MDUCT  - IMsoln.iI  a^d 

if  Msoln.i  > 0:  W “ IFR  = 1 and  ITO  = 2 

if  ^soln  i < 0:  IDir  ” -1  and  IFR  = 2 and  ITO  = 1 

if  Mg  q l u i ”0:  Iqir  “ 0,  there  is  no  forced  vent 

flow,  and  no  adjustments  are  made  to 
the  environments  in  rooms  1 and  2 


ii . 

if  nsoln 

- 2, 

then  IDIR  - IDIROld  and 

if 

Idir 

“ ^soln , i/I^soln , i 1 > then  MDUCT 

“ I^SOLN , 1 

if 

i D I R 

“ ol n ,2/1  ol n . 2 1 ’ then  MD u c T 

“ 1 0 L N , 2 

else : 


if  IDIR  - 1,  then  IFR  - 1 and  ITO  - 2 
else  IFR  - 2 and  ITO  - 1 
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4. 


Determine  the  number  and  properties  of  the  flow  slabs  leaving  the  room 
IFR  and  the  properties  of  the  duct  flow. 


a-  if  ySYSTEM,IFR  + HIFr/2  * yLAYER,IFR’  then  the  duCt  fl°W  is 

entirely  from  the  lower  layer.  Set  NSLABF  - 1.  Then: 

^SLABF ,1  “ ^DUCT 


^SLABF , 

, 02 

, 1 

‘ ^DUCT 

, 02 

“ MD  U C T * CL  , 

, 02  , 

IFR 

^SLABF , 

, K , 

1 “ 

U C T , 

K “ 

UCT  * CL , K , 

, IFR 

, K — 2 to  NpR0D 

Qslabf , 

, 1 

“ Qd 

UC  T “ 

UCT  ‘CpTL  , IFR 

fsLABF , 

, 1 

- Td 

UCT  “ 

tl. 

IFR 

ySLABF 

, 1 

- ys 

YSTEM , 

IFR 

b-  if  ySYSTEM.IFR  ’ HIFr/2  ~ yLAYER.IFR’  tben  tbe  duCt  flow  is 
entirely  from  the  upper  layer.  Set  NSLABF  - 1.  Then: 


^LABF  .1  = ^DUCT 

PSLABF,02,1  ” U C T , 0 2 “ UC  T * CU  , 0 2 , I FR 

PSLABF,K,1  “ PDUCT,K  “ MD U C T * CU , K . I FR » K “ 2 tO  NpRQD 

QsLABF.l  “ QdUCT  “ ^DUCT “Cp^U . IFR 

TSLABF,1  “ TDUCT  “ TU,IFR 

ySLABF.l  “ ySYSTEM.IFR- 

c.  else  the  duct  flow  is  from  both  the  upper  and  lower  layers.  Set 
NSLabf  " 2 and  iet  slab  1 come  from  the  upper  layer  and  slab  2 
come  from  the  lower  layer.  Compute  Tduct  from  Eq.  (8);  Mslabl.ifr 
from  Eq.  (9);  and  MgLABL  IFR  from  Eq.  (10).  Then: 

^SLABF ,1  “ MU , IFR 


^SLABF ,02,1 
^SLABF , K , 1 

Qslabf.i  “ 

NSLABF , 1 “ 
ySLABF , 1 “ 
L A B F , 2 ” 


“ MU . IFR ’CU , 02 , IFR 
“ MU , IFR  *CU , K , IFR , 
MU  , IFR  -CpTu , IFR 
TU  , IFR 

ySYSTEM.IFR  +bU,IFR 
ML  . I F R 


K - 


/2 


2 to  nprod 
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FSLABF,02,2  ML,IFR'CL,02,IFR 

FSLABF,K,1  = ML , IFR ’ CL , K , IFR , K = 2 to  NpR0D 


Qs  L A B F 

,2  = ML  , I F R ' Cp  TL  , 

IFR 

tslabf 

,2  = TL  , I F R 

ySLABF 

,2  = ySYSTEM , IFR 

, I F R 

/2 

FDUCT  , 

02  = FS  L A B F , 0 2 , 1 

+ fslabf 

,02,2 

P 

rD  U C T , 

K = FSLABF  , K , 1 + 

FSLABF  , K 

, 2 > K = 2 to  nprod 

Qduct 

QsLABF.l  + QsLABF,2 

tduct 

TD  U C T , IFR 

Determine  the  number  and  properties  of  the  flow  slabs  entering  room  ITO. 

a-  if  ^SYSTEM , IFR  + ^IFr/^  ~ ^LAYER , I FR  °r  ^SYSTEM , IFR  " ^IFr/^  — 

Ylayer  ifr  ’ t^ien  the  duct  flow  enters  entirely  the  upper  or  lower 
layer.’  Set  NSLABT  = 1.  Then: 

^slabt , 1 = mduct 

FSLABT,02,1  = PDUCT,02 

FSLABT,K,1  = FDUCT,K'  K = 2 tO  NpR0D 
Qslabt.i  = Qduct 

TSLABT  ,1  = f D U C T 
ySLABT.l  = ySYSTEM.ITO. 

b.  else  the  duct  flow  is  split  between  the  upper  and  lower  layer. 

Set  Nslabf  = 2.  Let  slab  1 enter  the  upper  layer  and  slab  2 enter 
the  lower  layer.  Let  Fy  and  FL  be  the  fraction  of  the  duct  flow 
going  to  the  upper  and  lower  layer,  respectively.  Then: 

Fu  = , i T o /Hi  T o 

Fl  = 1 - Fu 

ySLABT.l  = ySYSTEM.ITO  + N.ITO  /2 
ySLABT,2  = ySYSTEM.ITO  ' ^L.ITo/^ 

^SLABT,  1 = FU  'MDUCT 
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Qs  L A B T , 1 ^U  ’ Qd  U C T 


T, 


SLABT , 1 


Tr 


DUCT 


M = F .M 

SLABT,  2 rU  “DUCT 


Qs  L A B T , 2 ^U  ' Qd  U C T 


T, 


SLABT , 2 


= tt 


DUCT 


SUBROUTINES  USED 

FANRES  [2] 
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SUBROUTINE  VARIABLES 

All  nomenclature  in  the  subroutine  is  identical  to  the  nomenclature  used 
above  except  for: 

CL  , K , 1 > CL  , K , 2 ’ K “ 2 to  Nprod 


CL(K, 1) , CL(K, 2) ; K = 2 to  NPROD 


CU(K, 1) , CU(K, 2) ; K - 2 to  NPROD 
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CL , 0 2 . 1 • 

CL  , 0 2 , 2 

- 

CL(1 , 1) , CL(1 , 2) 

CU  , 02 , 1 » 

CU , 02 , 2 

- 

CU(1,1) , CU ( 1 , 2 ) 

CP 

- 

CP  [W- s/ (kg*K) ] 

g 

- 

G [m/s2 ] 

IR 

- 

DIRS12 

mduct 

- 

MDUCT  [kg/s] 

LAB  F >N  > 

^SLABT  >N 

- 

MSLBF(N),  MSLBT(N)  [kg/s] 

^DUCT 

- 

MFLOW  [kg/s] 

NSLABF  > N 

fSLABT 

- 

NSLBF,  NSLBT 

PD  U C T , 0 2 

- 

PDUCT(l)  [(kg  of  oxygen) /s] 

PD  U C T , K 

- 

PDUCT(K)  [(unit  of  product  K)/s] 

^SLABF , K , 

N ' ^SLABT  ’ K , N 

- 

PSLBF(N.K) ,PSLBT(N,K) 

[(unit  of  product  K)/s] 

S L AB  F , 02 , H ’ ^SLABT,02,N 

- 

PSLBF(N.l) ,PSLBT(N,1)  [(kg  of  oxygen)/s] 

Qduct 

QDUCT  [W] 

QsLABF , N > 

QsLABT  , N 

- 

QSLBF(N) , QSLBT(N)  [W] 

tduct 

- 

TDUCT  [K] 

L AB  F , N > 

^SLABT  , N 

- 

TSLBF(N) , TSLBT(N)  [K] 

TU  ’ 1 ’ V 

2 

- 

TU(1) , TU(2)  [K] 

tl>i.  Tl. 

2 

- 

TL(1) , TL(2)  [K] 

V* 

- 

V(J)  [ m3 /s ] 

^FLOOR , 1 > 

yFLOOR , 2 

- 

YFLOR(l) , YFLOR ( 2 ) [m] 

yLAYER , 1 » 

yLAYER , 2 

- 

YLAY(l) , Y1AY(2)  [m] 

ySLABF , N > 

ySLABT  , N 

- 

YSLBF(N) , YSLBT(N)  [m] 

ySYSTEM , 1 

> ySYSTEM , 2 

- 

YSYS(l) , YSYS(2)  [m] 
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’f  , J 

- 

DP(J)  [Pa  = kg/(m-s2)] 

’floor  , 1 > 

^Pfloor , 2 

PFLOR(l) , PFL0R(2)  [ Pa  - kg/(m • s2 ) ] 

’system  , 1 

> ^PsYSTEM,2 

PS (1) , PS (2)  [Pa  - kg/(m-s2)] 

, 1 » Pi.  ,2 

- 

DENL(l) , DENL(2)  [kg/m3] 

. 1 > P\]  , 2 

- 

DENU(l) , DENU(2)  [kg/m3] 

PREPARED  BY 

Leonard  Y.  Cooper  and  Glenn  P,  Forney 
January  1989 
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Figure  1.  The  fire -generated  environment  in  two  rooms  connected  by  a 
fan/ duct  system. 
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Figure  2 . Flow  entering  the  duct  in  room  IFR  when  the  layer  interface  there 
is  between  the  top  and  bottom  of  the  vent  opening. 
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Figure  3 . 


A sketch  of  a generic  fan- operating  curve  which  increases 
monotonically  in  volume -flow -rate  with  decreasing  cross-fan 
pressure  difference. 
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SUBROUTINE 


SUBROUTINE  VENTF( 


I 

ACC, 

FANPRS , 

FANVOL, 

JTERM,  CONL, 

CONU, 

I 

CP, 

NPROD , 

MXPRD, 

PFLOR,  RESIS, 

TL, 

I 

YLAY, 

YSYS, 

RHOL, 

RHOU,  HFVNT, 

LSTSOL 

0 

NSLBF, 

NSLBT , 

DIRS12 , 

0 

XMSLBF, 

XMSLBT, 

XMDUCT, 

QSLBF,  QSLBT, 

QDUCT , 

0 

PSLBF, 

PSLBT, 

PDUCT , 

TSLBF,  TSLBT, 

TDUCT , 

0 

f 

YSLBF, 

YSLBT 

Oc 

IMPLICIT 

DOUBLE 

) 

PRECISION  (A-H , O-Z) 

C*BEG 

C 


TU,  YFLOR , 


C***  VENTF 
C 
C 

C***  SUBROUTINE  ARGUMENTS 
C 


GENERIC  - THIS  ROUTINE  CALCULATES  THE  FLOWS  (MASS,  ENTHALPY 
AND  PRODUCT  OF  COMBUSTION)  THROUGH  A FAN  DUCT  SYSTEM. 


VOLUME  FLOWS  [M**3/S]  THAT 


INPUT 

ACC  - ACCURACY  OF  CALCULATION 
FANPRS , FANVOL  - ARRAYS  OF  PRESSURES  [P] 

DEFINE  THE  FAN  CURVE 

JTERM  - NUMBER  OF  TERMS  IN  THE  FANCURVE 
CONL  - CONCENTRATION  OF  ALL  PRODUCTS  FOR  LOWERLAYERS 

CONU  - CONCENTRATION  OF  ALL  PRODUCTS  FOR  UPPER  LAYERS 

CP  - SPECFIC  HEAT  AT  CONSTANT  PRESSURE  [W*S/(KG*K)] 

NPROD  - ACTUAL  NUMBER  OF  PRODUCTS 
MXPRD  - MAXIMUM  NUMBER  OF  PRODUCTS 
PFLOR  - PRESSURE  AT  FLOOR  [P] 

RESIS  - RESISTANCE  OF  FAN  SYSTEM 

TL.TU  - TEMPERATURE  OF  LOWER/UPPER  LAYERS  [K] 

YFLOR  - HEIGHT  OF  FLOORS  ABOVE  DATUM  ELEVATION  [M] 

YLAY  - HEIGHT  OF  LAYERS  ABOVE  FLOOR  [M] 

YSYS  - HEIGHT  OF  FAN  SYSTEM  ABOVE  FLOOR  [M] 

RHOL  - DENSITY  OF  LOWER  LAYERS  [KG/M**3] 

RHOU  - DENSITY  OF  UPPER  LAYERS  [KG/M**3] 

HFVNT  - WIDTH  OF  INLET/OUTLET  OF  FAN  SYSTEM  [M] 

LSTSOL  - LAST  SOLUTION  FOR  THIS  FAN  SYSTEM;  LSTSOL-1  —> 

PREVIOUS  SOLUTION  WAS  IN  FORWARD  DIRECTION,  LSTSOL  - 2 
— > PREVIOUS  SOLUTION  WAS  IN  REVERSE  DIRECTION 

OUTPUT 

NSLBF.NSLBT  - NUMBER  OF  SLABS  (1  OR  2)  IN  INLET/OUTLET 
DIRS 12  - DIRECTION  OF  FLOW  FOR  EACH  SLAB 
XMSLBF , XMSLBT , XMDUCT  - MASS  FLOW  FOR  FROM/TO,  BOTH  SLABS 
QSLBF , QSLBT , QDUCT  - ENTHALPY  FLOW  FOR  FROM/TO,  BOTH  SLABS 
PSLBF , PSLBT , PDUCT  - PRODUCT  FLOW  FOR  FROM/TO,  BOTH  SLABS 
TSLBF , TSLBT , TDUCT  - TEMPERATURE  FOR  FROM/TO,  BOTH  SLABS  [K] 
YSLBF, YSLBT  - HEIGHT  OF  FROM/TO  SLAB 
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C*END 


DIMENSION  FANPRS(*) , FANVOL(*) , CONL(MXPRD , 2) , CONU(MXPRD , 2) 
DIMENSION  PFLOR(2) , TL(2) , TU(2) 

DIMENSION  YFLOR(2) , YLAY ( 2 ) , YSYS(2),  RHOL(2) , RHOU(2) 
DIMENSION  PS (2) , HFVNT(2) 

DIMENSION  VFANMS (2) 

DIMENSION  XMSLBF(2) , XMSLBT(2) 

DIMENSION  QSLBF(2) , QSLBT(2) 

DIMENSION  PSLBF(2,MXPRD) , PSLBT(2 ,MXPRD) , PDUCT(MXPRD) 
DIMENSION  TSLBF(2) , TSLBT(2) , YSLBF(2) , YSLBT(2) 

DIMENSION  YVBOT(2) , YVTOP(2) 

DIMENSION  RES IS (2) 

INTEGER  DIRS12 
C 

DIMENSION  RHODCT (2) , TDCT(2) , FL(2) , FU(2) , HL(2) , HU(2) 
PARAMETER  (G-9.80D0) 

PARAMETER  (R=400 . 0D0/1 . 40D0) 

C 

DO  100  IROOM  -1,2 

YVTOP(IROOM)  - YSYS( IROOM)  + HFVNT ( IROOM) /2 . 0D0 
YVBOT( IROOM)  = YSYS( IROOM)  - HFVNT ( IROOM) /2 . 0D0 
C 

C***  CALCULATE  PS, I 
C 

IF ( YSYS ( IROOM) . GT . YLAY ( IROOM) ) THEN 
PS (IROOM)  - PFLOR( IROOM) 

$ -RHOL( IROOM) *G* (YLAY( IROOM) - YFLOR ( IROOM) ) 

$ -RHOU ( IROOM) *G* (YSYS ( IROOM) -YLAY( IROOM)) 

ELSE 

PS (IROOM)  - PFLOR( IROOM) 

$ -RHOL( IROOM) *G* (YSYS (IROOM) -YFLOR(IROOM) ) 

END  IF 
C 

C***  COMPUTE  RHODCT 
C 

IF ( YVTOP ( IROOM) . LE . YLAY( IROOM) ) THEN 
FL( IROOM)  = 1.0D0 
FU( IROOM)  - 0.0D0 

ELSE  IF (YVBOT( IROOM) .GE. YLAY (IROOM) ) THEN 
FL( IROOM)  - 0.0D0 
FU( IROOM)  - 1.0D0 
ELSE 

HU (IROOM)  - YVTOP (IROOM)  - YLAY (IROOM) 

HL( IROOM)  - YLAY( IROOM)  - YVBOT( IROOM) 

A - HL ( IROOM) *SQRT (RHOU (IROOM)) 

B - HU ( IROOM) *SQRT(RHOL( IROOM)) 

FL( IROOM)  - A/(A+B) 

FU( IROOM)  - B/(A+B) 

END  IF 

TDCT( IROOM)  - FL( IROOM) *TL( IROOM)  + FU( IROOM) *TU( IROOM) 
RHODCT (IROOM)  - RHOU ( IROOM) *TU( IROOM) /TDCT( IROOM) 
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100  CONTINUE 


C 

C***  CALCULATE  MASS  GOING  THROUGH  FAN 
C 

IF(JTERM. EQ. 1)THEN 

IF(FANVOL(l) . GE.0)THEN 

XMDUCT  - RHODCT ( 1 ) *FANVOL( 1 ) 

ELSE 

XMDUCT  - RHODCT (2) *FANVOL(l) 

END  IF 
ELSE 

ACC  -=  1.0D-8 

DPS  - PS(1)  - PS (2) 

DYS  - YSYS(l)  - YSYS(2) 

CALL  FANRES ( 

I ACC , FANPRS , FANVOL , JTERM , DPS , RES I S , DYS , RHODCT , 

0 VFANMS ,NSOLN 

+ ) 

IF(NSOLN. EQ. 1)THEN 
XMDUCT  = VFANMS (1) 

ELSE 

XMDUCT  - VFANMS (LSTSOL) 

CALL  MSGPRT( 'MULTIPLE  SOLUTIONS  IN  FANRES ',0) 

END  IF 
ENDIF 

C 

C***  DETERMINE  DIRECTION  OF  FLOW,  (THE  FROM  ROOM) 

C 

IF (XMDUCT . GT . 0 . 0D0 ) THEN 
IFROM  - 1 
DIRS12  - 1 

ELSE  IF (XMDUCT. LT.0.0D0) THEN 
IFROM  - 2 
DIRS12  - -1 
ELSE 

IFROM  - 1 
DIRS12  - 0 
RETURN 
ENDIF 

ITO  - 3 - IFROM 
XMDUCT  - ABS (XMDUCT) 

C 

C***  DETERMINE  THE  NUMBER  AND  PROPERTIES  OF  THE  FLOW  SLABS  LEAVING 
C THE  ROOM  IFROM  AND  THE  PROPERTIES  OF  THE  DUCT  FLOW 
C 

IF(YVTOP( IFROM) . LE . YLAY ( IFROM) ) THEN 
NSLBF  - 1 

XMSLBF(l)  - XMDUCT 

QDUCT  - XMDUCT*CP*TL( IFROM) 

QSLBF(l)  - QDUCT 
TDUCT  - TL( IFROM) 

TSLBF(l)  - TDUCT 
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YSLBF(l)  - YSYS(IFROM) 

DO  210  I PROD  - 1,  NPROD 

PDUCT(IPROD)  - XMDUCT*CONL(IPROD , I FROM) 

PSLBF(1, IPROD)  = PDUCT(IPROD) 

210  CONTINUE 

ELSE  IF(YVBOT(IFROM) . GE . YLAY ( IFROM) ) THEN 
NSLBF  = 1 

XMSLBF(l)  - XMDUCT 

QDUCT  - XMDUCT*CP*TU( IFROM) 

QSLBF(l)  = QDUCT 
TDUCT  - TU( IFROM) 

TSLBF(l)  - TDUCT 
YSLBF(l)  - YSYS( IFROM) 

DO  220  IPROD  - 1,  NPROD 

PDUCT(IPROD)  - XMDUCT*CONU( I PROD, IFROM) 

PSLBF(1, IPROD)  = PDUCT(IPROD) 

220  CONTINUE 
ELSE 

NSLBF  - 2 

XMLIFR  = XMDUCT*FL( IFROM) 

XMUIFR  - XMDUCT*FU( IFROM) 

TDUCT  - TDCT( IFROM) 

XMSLBF(l)  - XMUIFR 

QSLBF(l)  - XMUIFR*CP*TU( IFROM) 

TSLBF(l)  - TU( IFROM) 

YSLBF(l)  - YLAY (IFROM)  + HU ( IFROM) /2 . 0D0 
C 

XMSLBF(2)  - XMLIFR 

QSLBF(2)  - XMLIFR*CP*TL( IFROM) 

TSLBF(2)  - TL( IFROM) 

YSLBF(2)  = YLAY( IFROM)  - HL( IFROM) /2 . 0D0 
QDUCT  - QSLBF(l)  + QSLBF(2) 

DO  250  IPROD  - 1,  NPROD 

PSLBF(1, IPROD)  - XMUIFR*CONU( I PROD, IFROM) 

PSLBF(2 , IPROD)  - XMLIFR*CONL( IPROD , IFROM) 

PDUCT(IPROD)  - PSLBF(1, IPROD)  + PSLBF(2 , IPROD) 

250  CONTINUE 

END  IF 
C 

C***  DETERMINE  THE  NUMBER  AND  PROPERTIES  OF  THE  FLOW  SLABS  ENTERING 
C THE  ROOM  ITO  - 3 - IFROM  AND  THE  PROPERTIES  OF  THE  DUCT  FLOW 
C 

ITO  - 3 - IFROM 

IF(YVTOP(ITO) . LE. YLAY (ITO) . OR. YVBOT(ITO) . GE. YLAY ( ITO )) THEN 
NSLBT  - 1 

XMSLBT(l)  - XMDUCT 
QSLBT(l)  - QDUCT 
TSLBT(l)  - TDUCT 
YSLBT(l)  - YSYS(ITO) 

DO  310  IPROD  - 1,  NPROD 

PSLBT(1, IPROD)  - PDUCT(IPROD) 


VENTF 


23 


310  CONTINUE 
ELSE 

NSLBT  - 2 

FFU  - HU(ITO) /HFVNT(ITO) 

FFL  - l.ODO  - FFU 

YSLBT(l)  - YLAY(ITO)  + HU(IT0)/2 . ODO 

YSLBT(2)  - YLAY(ITO)  - HL(IT0)/2 . ODO 

XMSLBT(l)  = FFU*XMDUCT 

XMSLBT(2)  - FFL*XMDUCT 

QSLBT(l)  - FFU*QDUCT 

QSLBT(2)  - FFL*QDUCT 

TSLBT(l)  - TDUCT 

TSLBT(2)  - TDUCT 

DO  320  IPROD  - 1,  NPROD 

PSLBT(1, IPROD)  - FFU*PDUCT ( I PROD ) 
PSLBT(2, IPROD)  - FFL*PDUCT( IPROD) 
320  CONTINUE 
END  IF 
RETURN 
END 
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VENTHP 


Calculation  of  the  Flow  of  Mass,  Enthalpy,  Oxygen  and  Other 
Products  of  Combustion  Through  a Vertical,  Constant-Width  Vent  in 
a Wall  Segment  Common  to  Two  Rooms 


DESCRIPTION 


Consider  an  instant  of  time  during  the  simulation  of  a multi-room  compartment 
fire  environment.  This  algorithm  can  be  used  to  calculate  characteristics  of 
the  flow  of  mass,  enthalpy,  oxygen,  and  other  products  of  combustion  through  a 
constant -width  vent  located  in  a wall  segment  common  to  any  two  rooms  of  the 
compartment  [1] . 

Designate  the  two  rooms  of  interest  as  rooms  1 and  2 and  refer  to  Figure  1 
which  depicts  the  room  fire  environments.  It  is  assumed  that  the  temperature, 
density,  concentration  of  oxygen  and  of  other  products  of  combustion  in  the 
upper  and  lower  layer  of  each  room  are  specified.  Also  specified  in  each  room 
are:  the  elevation  above  a datum  elevation  of  the  floor,  and  the  smoke -layer 
interface ; the  pressure  at  the  floor  of  each  room  above  a specified  datum 
pressure  and  the  absolute  pressure  there. 

In  general,  within  the  elevation  interval  of  the  common  wall  segment  the  room 
1-to-room  2 pressure  difference  can  be  described  by  at  most  three,  piecewise- 
continuous,  linear  functions  of  elevation.  These  will  have  endpoints  at  the 
bottom  and  top  of  the  common  wall  segment  and  at  the  elevations  of  the 
interfaces  which  separate  the  upper  and  lower  layers  of  the  two  rooms 
(provided  these  interfaces  also  lie  in  the  elevation  interval  of  the  wall 
segment).  This  is  depicted  in  Figure  2.  As  illustrated  in  this  figure,  it  is 
possible  for  each  of  the  pieces  of  the  piecewise  linear  functions  to  change 
sign  somewhere  within  its  operable  elevation  range.  The  elevations  (no  more 
than  four)  and  corresponding  pressure  differences  at  the  endpoints  of  these 
pieces  are  required  for  the  present  calculation.  Also  required  are  the 
elevations  (no  more  than  three  in  addition  to  the  above)  where  the  pieces 
change  sign  (i.e.,  elevations  where  the  room  1-to-room  2 pressure  difference 
is  zero) . 

It  is  assumed  that  the  number  and  values  of  all  of  the  above  elevations  and 
the  corresponding  absolute  pressures  in  each  room  as  well  as  the  pressure 
differences  have  been  determined  previously  and  that  these  are  available  as 
input  data  (e.g.,  from  use  of  C0MWL1  [2]).  NWELEV  < 7 is  the  number  of  the 
unique  elevations  and  yWELEV>N,  Pwi,n>  Pw2,n-  and  aPwi,w2,n  = Pwi.h  ' Pw2,n>  n 
= 1 to  NWELEV , are  the  values  of  the  elevations,  increasing  uniformly  with  N, 
the  corresponding  absolute  pressures  in  room  1 and  2,  and  the  pressure 
differences,  respectively. 

Consider  a vent  located  in  the  common  wall  segment.  The  type  of  vent 
considered  here  is  of  constant  width  and  has  a specified  top  and  bottom 
elevation,  yVTOp  and  yVB0T , respectively,  where  yVTOp  - yVB0T  > 0.  This 
algorithm  first  computes  the  pressures  and  room  1-to-room  2 pressure 
difference  at  yVT0P  and  yVB0T  and  forms  new  arrays  yN , px ( N , P2 , N > APi , 2 , n by 
sorting  the  results  into  the  original  input  arrays,  yWELEv  n>  Pwi  n>  Pw2,N’ 
and  Apwl  W2  N.  All  elements  of  the  new  yN  above  yVT0P  and  below  yVB0T , along 
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with  corresponding  px  N , P2 , n > and  APi  2 n are  ^en  eliminated.  The  result 
are  the  uniquely  determined  elevations  of  the  endpoints  and  zeros  of  the 
piecewise-continuous , linear,  pressure-difference  function  of  elevation 
bounded  by  yVT0P  and  yVB0T • NVELEV  < 7,  the  total  number  of  these  endpoints 
and  zeros,  are  found.  The  result  is  given  as  the  output,  yN , px  N,  p2  N, 

Api  2 N,  N = 1 to  NVELEV.  These  are  the  elevations,  increasing  uniformly  with 
N,  and  corresponding  absolute  pressures  in  room  1 and  2 and  pressure 
differences,  respectively.  Note  that  yx  = yVB0T  and  ^nelev  = Yviop  • 

The  gases  in  all  layers  of  the  two  rooms  are  modeled  as  the  same  perfect  gas. 
Recommended  properties  are  those  of  air. 

Associated  with  each  yN , yN+1  pair  is  a "slab"  of  flow,  designated  as  slab  N, 
which  is  uniform  in  composition  and  which,  by  virtue  of  the  cross -vent 
pressure  difference,  is  driven  between  the  two  rooms.  The  number  of  slabs  of 
flow  passing  through  the  vent  is  NSLAB  = (NVELEV  - 1)  <7.  Flow  is  from  room 
1 to  room  2,  or  from  room  2 to  room  1 depending  on  whether  either  of  Apx  2 N 
or  Apx  2 N+1  is  positive  or  negative,  respectively  (one  of  the  pair  can  be 
zero,  but  it  is  impossible  for  them  to  be  of  different  signs).  The  direction 
of  the  flow  is  determined  and  identified  by  the  value  IDIR  N.  For  flow  from 
room  1 to  room  2 or  from  room  2 to  room  1,  the  value  of  IDIR  N is  1 or  -1, 
respectively.  IDIR  n is  0 if  the  room-to-room  flow  is  identically  zero  (this 
can  only  occur  if  both  Apx  2 N and  Apx  2 N+1  are  identically  zero). 

Depicted  in  Figure  3 are  possible  distributions  of  cross-vent  pressure  for  a 
single  vent  and  the  direction  of  the  velocity  in  each  of  its  one -to -six 
uniform-property  flow  slabs.  The  details  of  such  a plot  for  a particular  vent 
at  a given  instant  of  time  is  determined  from  the  arrays  yN , Apx  2 N , and 

^DIR , N • 

For  each  flow  slab,  N,  the  algorithm  finds  the  density,  ps LAb,n>  total  mass 
flow  rate,  MgLAB  N,  absolute  temperature,  TSLAB  N,  enthalpy  flow  rate, 

Qslab.n-  oxygen  concentration,  cSLAB>02  N,  oxygen  mass  flow  rate,  P02(slab.n> 
and  the  concentrations  and  flow  rates  of  the  rest  of  the  NPROD  products  of 
combustion  being  tracked  by  the  simulation,  cSLAB  K N and  PK  SLAB  N , 
respectively,  K = 2 to  NPROD. 

Ts lab  n and  lab  n are  t^e  Properties  of  the  slab  as  it  enters  the  receiving 
room.  These  are  computed  with  the  assumption  that  the  flow  process,  from 
stagnation  in  the  source  room  to  the  point  where  the  flow  has  decelerated  in 
the  receiving  room  (without  significant  mixing),  is  one  of  constant  enthalpy. 
Thus,  Tslab  n is  the  absolute  temperature  of  the  slab's  source  layer. 

Also  determined  by  the  algorithm  are  the  elevations  above  the  datum,  ySLABN> 
of  the  centroids  of  the  distributions  of  the  flow  slabs'  entering  horizontal- 
velocity,  or  -momentum  profiles,  i.e., 


ryN+1  -^n+i  yN+i  ryN+i 

ySLAB,N  3 J pw ydy  / J ^vdy  * j ^y  / J vdy  > N = 1 to  nslab 
yN  yN  yN  yN 


(i) 
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where  y is  elevation  above  a common  datum,  and  p and  v = v(y)  are  the  density 
and  velocity  of  the  flow  as  it  enters  the  receiving  room.  In  Eq . (1),  the 
variation  of  p with  velocity  is  neglected  and  p is  approximated  as  being 
uniform  through  the  slab.  To  the  extent  that  the  horizontal  penetration  of  a 
slab  into  the  receiving  room  is  dependent  on  its  horizontal  momentum,  ySLAB  N 
may  be  preferable  to  other  characteristic  elevations  as  a measure  of  the 
elevation  of  its  flow  as  it  is  projected  into  the  receiving  room  by  its 
horizontal  momentum. 

The  present  algorithm  is  applicable  for  arbitrary  values  of  room  pressures  and 
cross-vent  pressure  differences,  i.e.,  whether  or  not  compressibility  plays  a 
significant  role  in  the  vent  flow.  The  flow  calculation  used  here  is  from 
[3].  This  is  an  extension  to  the  incompressible  vent  flow  calculation 
presented,  for  example,  in  [4]. 


OUTPUT 


CS  L A B , 0 2 , N ’ N 1 tO  Nslab 

Concentration  of  oxygen  in  slab  N between  the 
elevations  yN  and  yN+1.  [ (leg  of  oxygen)/(kg  of  slab 

flow) ] 

cs L ab , K , n > N - 1 to  NSLAB;  K - 2 to  NPROD 

Concentration  of  product  of  combustion  K in  slab  N 
between  the  elevations  yN  and  yN  + 1.  [(unit  of 
product)/(kg  of  slab  flow)] 


1 to  Nc 


A measure  of  the  direction  of  the  room  1- to -room  2 
flow  in  slab  N,  between  yN  and  yN+1 . For  flow  from 
room  1 to  room  2 or  from  room  2 to  room  1,  the  value 
of  IDIR>N  is  1 or  -1,  respectively.  IDIRN  is  0 if 
the  room-to-room  flow  is  identically  zero. 


APi , 2 . n ; N - 1 to  Nvelev 

Pi  n ‘ P2  n>  i-e->  pressure  in  room  1 - pressure  in 
room  2 at’ elevation  yN . If  Apx  2 N >0,  - 0,  or  < 0 
then  flow  through  a vent  at  this  elevation  will  be 
from  room  1 to  room  2,  zero,  or  from  room  2 to  room  1, 
respectively.  [N/m2  = kg/(m-s2)  = pascal] 
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; N = 1 to  N, 


M 


SLAB  , N 


N 


SLAB 


SLAB 


Total  mass  flow  rate  in  slab  N between  the  elevations 
yN  and  yN+1 . [kg/s] 


Number  of  slabs  of  flow,  < 6,  between  the  bottom  and 
top  of  the  vent. 


^VELEV 


Number  of  unique  elevations,  < 7,  which  identify:  1) 
endpoints  of  the  continuous,  linear  pieces  of  the  room 
1-to-room  2 pressure  difference  vs  elevation;  and/or 
2)  elevations  of  zero  room  1-to-room  2 pressure 
difference,  where  all  these  elevations  are  at  or 
between  the  top  and  bottom  of  the  vent. 


P02 , SLAB  , N > 


N - 1 to  Nslab 


Mass  flow  rate  of  oxygen  in  slab  N between  the 
elevations  yN  and  yN+1.  [(kg  of  oxygen)/s] 


^k,slab,n>  ^ 1 to  NSLAB ; K - 2 to  NpR0D 

Flow  rate  of  product  of  combustion  K in  slab  N 
between  the  elevations  yN  and  yN+1.  [(unit  of  product 
K)/s  ] 


Pi,N  [p2iN] ; N = l to  nvelev 

Absolute  hydrostatic  pressure  in  room  1 [room  2]  at 
elevation  yN . [Pa  = kg/(m-s2)] 


Qslab,n:  N - 1 to  Nslab 

Total  enthalpy  flow  rate  in  slab  N based  on  absolute 
temperature  between  the  elevations  yN  and  yN+1 . [W] 


^slab.N’  N - 1 to  Nslab 

Density  of  the  flow  in  slab  N between  the  elevations 
yN  and  yN+1  as  it  enters  the  receiving  room  [kg/m3]. 
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; N = 1 to  N, 


^SLAB.N 


SLAB 

Absolute  temperature  of  the  flow  in  slab  N between  the 
elevations  yN  and  yN+1  as  it  enters  the  receiving 
room.  [K] 


ySLAB,N;  N '1  t0  NSLAB  ’ anC*  ySLAB.N  < ySLAB  , N+l 

Elevations  above  the  datum  of  the  centroids  of  the 
distribution  of  slab  N's  entering  horizontal- 
velocity,  or  -momentum  profiles  as  determined  by  Eq . 
(1) . These  are  useful  as  a single  measure  of  the 
elevation  of  slab  N as  it  is  projected  into  the 
receiving  room,  [m] 


yN;  N = 1 to  Nvelev;  yN  < yN+1 

The  elevations  above  the  datum  elevation  which 
identify  the  bounds  of  the  NVELEV  uniform-composition 
slabs  of  flow  through  a vent,  [m] 


INPUT 


Av 


Area  of  the  vent  [m2 ] . 


cl , k , i tcu , k , i 1 > K - 2 to  NpR0D ; I 1 or  2 

Concentration  of  product  of  combustion  K in  lower 
[upper]  layer  of  room  I.  [(unit  of  product  k)/(kg  of 
layer) ] 


CL,02,I  tCU,02,l]’  ^ — 1 or  2 

Concentration  of  oxygen  in  lower  [upper]  layer  of  room 
I.  [(kg  of  oxygen)/(kg  of  layer)] 


Specific  heat  at  constant  pressure  of  the  vent  flow 
[W-s/(kg-K)]  (suggest  103  W-s/(kg-K)  for  air  as 
default) . 
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aPwi,w2,n ; N = 1 to  Nwelev 

Pwi  n ' Pw2  n ’ i-e-.  pressure  in  room  1 - pressure  in 
room  2 at  elevation  yWELEVjN.  If  aPwi.w2,n  > °-  = °> 
or  < 0 then  flow  through  a vent  at  this  elevation  will 
be  from  room  1 to  room  2,  zero,  or  from  room  2 to  room 
1,  respectively.  [Pa  = kg/(m-s2)] 

Nwelev 

Number  of  unique  elevations,  < 7,  which  identify:  1) 
endpoints  of  the  continuous,  linear  pieces  of  the  room 
1-to-room  2 pressure  difference  vs  elevation;  and/or 
2)  elevations  of  zero  room  1-to-room  2 pressure 
difference,  where  all  these  elevations  are  at  or 
between  the  horizontal  extremes  of  the  wall  segment 
common  to  rooms  1 and  2 . 


N 

1 r R 0 D 


Number  of  products  of  combustion,  including  oxygen, 
being  tracked  in  the  simulation. 


Pwi.N  [Pw2.r] ; N = 1 to  nwelev 

Absolute  hydrostatic  pressure  in  room  1 [room  2]  at 
elevation  yWEEEV , n • tPa  = kg/(m-s2)] 


Pdatum 


Datum  absolute  pressure.  [Pa  - kg/(m-s2)] 

5Pfloor,i ! I = 1 or  2 

Pressure  at  floor  in  room  1 [room  2]  above  the  datum 
pressure.  [Pa  = kg/(m-s2)] 


Vi  [Tt>1];  I = 1 or  2 

Absolute  temperature  of  the  upper  [lower]  layer  in 
room  I.  [K] 


ywELEV.N*  N - 1 to  NWELEV  ; and  yWELEV,N  yWELEV,N  + l 

Values  of  the  NWELEV  elevations  above  the  datum  eleva- 
tion. [m] 
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^FLOOR , I ^LAYER , I 


] ; I = 1 or  2 


Elevation  of  floor  [layer]  in  room  I above  the  datum 
elevation,  [m] 


yVTOP  [yVBOT ] 


Elevation  of  the  top  [bottom]  of  the  vent  above  the 
datum  elevation,  [m] 


Pl.i  [Pv, 1 1 ; 1 - 1 or  2 

Density  of  lower  [upper]  layer  in  room  I.  [kg/m3] 


Error  tolerance  for  5pFLOor,i-  If  Perror.i  is  the 
uncertainty  in  SpFL00RI,  I - 1 or  2 , then  pERR0R>I 
satisfies 

IPerror.iI  < €pPo  + I 5Pfloor  ’I  I eP 

where  p0  *=  1.0  Pa.  The  first  term  is  based  on  an 
absolute  error  tolerance  and  dominates  the  above  error 
bound  when  |SpFL00R  j|  is  less  than  1.0  Pa.  The 
second  term  is  a relative  error  tolerance  and 
dominates  when  |5pFL00R  j|  is  greater  than  1.0  Pa. 

The  error  tolerance,  ep , is  chosen  to  be  consistent 
with  the  tolerance  specified  for  the  pressure  equation 
in  an  ODE  solver  that  calls  this  routine. 

CALCULATION 

1.  Place  yWELEV,N>  Pwi,N>  Pw2,N>  APwi(W2.N  int0  yN  . Pi  , N > ?2  , N > APl , 2 , N > 
eliminating  all  array  sets  with  yN  below  yVB0T  and  above  yVtop • 

2.  Add  yVB0T  and  yVTOp  to  yN . 

3.  Using  the  algorithm/subroutine  DELP  [5],  calculate  the  absolute 
pressures  in  the  rooms  1 and  2 and  their  pressure  differences  at  the 
elevations  yVB0T  and  yVT0P . Store  the  results  at  appropriate  locations 
in  the  arrays  Pl  N , p2  N , Apx f 2 N . 

4.  Using  the  algorithms/subroutines  HSORT  [6,7]  and  RMVDUP  [8],  sort  the 
height/pressure -pair/pressure -difference  list  and  remove  any  entries 
with  duplicate  elevations. 
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5.  Set  the  number  of  slabs,  NSLAB , to  one  less  than  the  number  of 
elevations  in  yN . 

6.  Determine: 


IDIR  N - -1,  0,  1 according  to  whether 

(2) 

sign[Apx  2 N + Api ( 2 , n+ l ] < 0*  “ 0,  or  > 0,  respectively. 

7.  Determine: 


Tslab  n - absolute  temperature  of  slab  N's  source  layer  (3) 


PSlab  n “ (density  of  slab  N's  source  layer)" 

(absolute  pressure  at  slab  base  in  receiving  room)/ 
(absolute  pressure  at  slab  base  in  source  room) 

(4) 


cslab  02  n = concentration  of  oxygen  in  slab  N's  source  layer  (5) 


‘slab  k H “ concentration  of  product  of  combustion  K in  slab  N's 

source  layer,  K = 2 to  NpR0D 


(6) 


8.  For  nonzero- flow  slabs,  i.e.,  where  max(|Ap1  2 N|,|Ap1  2 N+1|)  > 0, 
set: 


p - 

LAB  . 

N > 

yT  - 

yN+i; 

yB  - 

yN; 

pt  - 

lAPl,2 

, N + 1 1 

; pb 

Pabst 

— max 

•[Pl.HH 

Pabsb 

- max 

[Pl.N 

• p2 

(7) 


where  the  subscripts  B and  T refer  to  the  bottom  and  top  of  the  slab, 
respectively.  Then  define 


VENTHP 


8 


e (Pt  + Pb )/ (Pabst  + Pa b s B ^ > x — 1 - e 
Following  [3]  compute  the  flow  coefficient  C(x)  from 
C (x)  = 0.85  - 0 . 17x  = 0.68  + 0.17e, 


(8) 


(9) 


and  w(x)  from 


w(x) 


f (x) 


1 - [3/(47) ]e  if  0 < e < 10' 5 
f (x) / [ 2e ] 1 7 2 if  1 > e > 10'5 

{ [27/(7-l) ]x2/7 [l-x(7'1)/7] }1/2  if  e < 1 - [2/(7+l) ]7/ (7_1) 
(7[2/(7+l)](7+1)/ (7-1) )1/2  if  e > 1 - [2/(7+l)]7/ (7_1) 


(10) 


(11) 


where  7,  the  ratio  of  specific  heats  of  the  vent  flow  gas,  is  taken  to 
be  that  of  air,  1.40. 


Define  p and  pCUT  by 

0 if  pT1/z  + pB1/2  = 0 

[PT  + (PbPt>1/2  + PbI/IPi1'2  + Pb 1 f 2 ] otherwise 

pHt  - [ep-MAX(1.0  pascal,|SpFL00R  J.ISppboob.bI)]1'2 

Compute  the  numerical  damping  factor,  FNOise>  from 

FN  OISE  = 1-°  ' eXP('P1/2/Pc(jT) 

a.  Determine  (from  Eq.  (17)  of  [3]): 

^SLAB.N  = FN  0ISE^^^W^X^  [(8p)^^/3]  Av 

t [yx  ■ yB]/[yvTop  _ yvBOT^’P1^2 


(12) 


(13) 


(14) 


(15) 
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Note  that  e -*  0 corresponds  to  the  incompressible  limit  and  that 
under  this  limit  the  mass  flow  rate  result  of  Eq.  (15)  is  identi- 
cal to  the  traditional  mass  flow-rate  result  of  Eq.  (33)  of  [4] 
(when,  as  recommended  in  [4] , the  vent  flow  coefficient  there  is 
taken  to  be  0.68). 

The  term,  FN0ISE , is  designed  to  damp  out  the  numerical  noise 
(error)  in  the  calculated  value  for  p that  would  otherwise  be 
dominant  in  Eq.  (15)  when  the  cross -vent  pressure  difference  terms 
pT , pB  are  small  relative  to  the  maximum  of  1.0  pascal  and  the 
calculated  floor  pressures,  5pFL00R1,  fipFL00Rt2.  The  term,  pCUT , 
is  an  estimate  of  how  small  the  maximum  of  (pT , pB ) must  be  to 
retain  a few  digits  of  accuracy  in  the  calculation  of  p.  When  the 
calculated  value  of  p is  smaller  than  pCUT , this  value  and, 
therefore,  the  value  of  MsLAB  N will  likely  contain  noise  which 
should  be  damped.  FN0ISE  is  constructed  to  tend  towards  1 when  pT 
or  pB  are  large  relative  to  pCUT  and  tends  towards  0 when  pT  and 
pB  are  small  relative  to  pCUT . 

As  an  example  of  the  importance  of  FN0ISE  consider  pB  for  a vent 
whose  bottom  is  at  the  floor  and  connects  two  rooms  with  equal 
floor  height.  Then  for  the  lowest  vent  flow  slab  pB  is  given  by 


Pfi  ” I^PFL00R,1  " ^PFL00R.2l  “ I 1,2,3.  I 

Both  5pFL00R  1 and  5pFL00R  2 are  quantities  that  are  typically  es- 
timated by  an  ODE  solver  to  only  a certain  user- specif ied 
accuracy,  say  ep  - 10' 6 . In  this  case,  if  the  first  6 digits  of 
5pFLooR  l an<*  ^Pfloor  2 match  then  the  difference,  pB  , will 
consist  entirely  of  random  noise  or  error.  Under  the  not- 
unreasonable  worst-case  scenario,  the  same  level  of  error  can  be 
expected  in  the  calculated  value  of  pT . The  numerical  errors 
present  in  pB  and  pT  leads  to  a similar  error  in  p which  is 
amplified  when  computing  MgLAB  N in  Eq.  (15).  This  is  amplified 
further  by  Cp*TSLABN  when  calculating  QSLAB>N  per  Eq.  (18)  below. 
The  term  QSLAB  N is  a component  of  the  differential  equation  for 
pressure.  So,  if  the  noise  in  Qslab  n dominates  the  pressure 
equation  the  ODE  solver  will  reduce  the  time  step  resulting  in 
solution  inefficiency  or  worse,  a solver  failure. 

With  MsLAB  n , now  calculate 

*02,  SLAB,  N “ CSLAB  , 02  , N *^SLAB  , N 

*K,SLAB,N  “ CSLAB , K , N '^SLAB , N > K - 2 tO  NpR0D  (17) 

Qslab, n " ^p ’^slab , n ‘^slab , n (18) 
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From  indicated  integrations  in  Eq . (1)  determine  (see  Section 
3. 3.4. 5 of  [1]): 


c . 


ySLAB.N  = [ (3/5)  (pT  5 / 2 - pB5/2)(yT  - yB )/ (pT3 /2  - Pb3/2) 

+ (pTyB  " pByT)]/(pT  ■ pb ) 


(19) 


In  order  to  increase  accuracy  when  calculating  the  above  bracketed 
term  when  pT  is  close  to  pB , the  subroutine  computes  ySLAB  N 
according  to: 


<yB  + yT)/2  if  !pb  - Pxl  * 0.001 

[ (3/5) (pT2  + Pt3/2Pb1/2  + PTPB 

(21 

+ pT1/2pB3/2+  pB2 ) (yT  - yB )/ 

(pT  + pt1/2pb1/2  + pb)  + (pTyB  - pByi)]/ 

(pT  - pB ) otherwise 


SUBROUTINES  USED 

DELP  [5] 

HSORT  [6,7] 
RMVDUP  [8] 
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SUBROUTINE  VARIABLES 


All  nomenclature  in  the  subroutine  is  identical  to  the  nomenclature  used  above 
except  for: 

Av  - AVENT  [m2] 

Cp  - CP  [ W • s/(kg • K) ] 

cs L a B , o 2 , N ■ CSLAB(l)  [(kg  oxygen)/ (kg  of  layer)] 

cs  lab , K , n ' CSLAB(K) 

[(unit  of  product  K)/(kg  of  layer)] 
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CL  , K , 1 ' CL  , K , 2 

' 

CONL(K.l),  CONL(K, 2) 

[(unit  of  product  K)/(kg  of  layer)] 

CU  , K , 1 » CU  , K , 2 

- 

CONU (K, 1) , CONU(K, 2) 

[(unit  of  product  K)/(kg  of  layer)] 

CL  , 02  , 1 ' CL , 02 , 2 

- 

C0NL(1 , 1) , CONL( 1 , 2) 

[(kg  of  oxygen)/ (kg  of  layer)] 

CU , 02 , 1 * CU , 02 , 2 

- 

C0NU(1,1),  C0NU(1,2) 

[(kg  of  oxygen)/(kg  of  layer)] 

f 

- 

FF  [dimensionless] 

^DIR  , N 

- 

DIRS12  [dimensionless] 

^SLAB  , N 

- 

XMSLAB  [kg/s] 

nslab  ’ nvelev  > ^WELEV 

- 

NSLAB,  NVELEV,  NMELEV  [dimensionless 

PK , SLAB  , N 

- 

PSLAB(K,N)  [(unit  of  product  K)/s] 

P02 , SLAB  , N 

- 

PSLAB(1,N)  [(kg  oxygen) /s] 

Pdatum 

- 

PDATM  [Pa  = kg/(m-s2)] 

Pw 1 , N > Pw  2 , N 

- 

PAB1 , PAB2  [Pa  = kg/(m-s2)] 

Pl,R»  P2,N 

- 

PVAB1 , PVAB2  [Pa  - kg/(m-s2)] 

Pabst  > Pabsb 

- 

PABST,  PABSB  [Pa  - kg/(m-s2)] 

Pb.  PT 

- 

PI , P2  [Pa  - kg/(m- s2 ) ] 

QsLAB  , N 

- 

QSLAB  [W] 

^SLAB  , N 

- 

TSLAB  [K] 

TU , 1 > TU  , 2 

- 

TU(1) , TU(2)  [K] 

TL,1.  TL,2 

- 

TL(1) , TL(2)  [K] 

w 

- 

W [dimensionless] 

X 

- 

X [dimensionless] 

^FLOOR  , 1 > yFLOOR , 2 

- 

YFLOR(l) , YFLOR ( 2 ) [m] 

Flayer , 1 » yLAYER , 2 

- 

YLAY(l) , YLAY(2)  [m] 

ySLAB  , N 

- 

YSLAB  [m] 
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yVTOP > yVBOT 

yWELEV 

yB.  yT 
yN 

7 

APW 1 , W 2 , N 
APl , 2 , N 

^Pfloor.i,  ^PfLOOR,2 
e 

S 

P 

P-L  , 1 - Ph  , 2 
P\J  , 1 ' P\]  , 2 

LAB  , N 

PREPARED  BY 

Leonard  Y.  Cooper  and  Glenn  P. 
March  1988 


YVTOP,  YVBOT  [m] 

YELEV  [dimensionless] 

Yl,  Y2  [m] 

YVELEV  [m] 

.1.40 

DP1M2  [Pa  = kg/(m-s2)] 

DPV1M2  [Pa  = kg/(m- s2 ) ] 

PFLOR(l) , PFL0R(2)  [Pa  - kg/(m- s2 ) ] 
EPS  [dimensionless] 

EPSP  [dimensionless] 

DEN  [kg/m3] 

DENL(l) , DENL(2)  [kg/m3] 

DENU(l) , DENU ( 2 ) [kg/m3] 

RSLAB  [kg/m3] 
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Figure 


1.  The  fire -generated  environment  in  two  adjacent  rooms. 
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Figure  2.  Possible  distributions  of  cross -wall  pressure  difference,  Apx  2 
for  a common- wall  segment. 
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Figure  3. 


Possible  distributions  of  cross-vent  pressure  difference,  Apj  2 
for  a single  vent  and  flow  directions  of  its  one-to-six  uniform 
property  flow  slabs. 
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SUBROUTINE 


SUBROUTINE  VENTHP ( 

I YFLOR , YLAY , 

I TU , TL , DENL , DENU , PFLOR , YVTOP , YVBOT , AVENT , YELEV , DP1M2 , NELEV , 

I CP , PDATM , PAB1 , PAB2 , CONL , CONU , 

I NPROD , MXPRD , MXSLAB.EPSP, 

OCSLAB , PSLAB , PV1 , PV2 , QSLAB , 

0 DIRS12 , DPV1M2 , RSLAB , TSLAB , YSLAB , YVELEV , XMSLAB , NVELEV , NSLAB) 

IMPLICIT  DOUBLE  PRECISION  (A-H.O-Z) 

C*BEG 

C***  VENTHP  GENERIC  - CALCULATION  OF  THE  FLOW  OF  MASS,  ENTHALPY,  OXYGEN 
C AND  OTHER  PRODUCTS  OF  COMBUSTION  THROUGH  A VERTICAL, 

C CONSTANT -WIDTH  VENT  IN  A WALL  SEGMENT  COMMON  TO  TWO  ROOMS. 

C THE  SUBROUTINE  USES  INPUT  DATA  DESCRIBING  THE  TWO-LAYER 

C ENVIRONMENT  IN  EACH  OF  THE  TWO  ROOMS  AND  OTHER  INPUT  DATA 

C CALCULATED  IN  SUBROUTINE  COMWL1 . 

C 

C***  SUBROUTINE  ARGUMENTS 
C 

C INPUT 

C 

C YFLOR 
C YLAY 
C TU 
C TL 
C DENL 
C DENU 
C PFLOR 
C 

C YVTOP 
C YVBOT 
C AVENT 
C YELEV 
C DP1M2 
C 

C NELEV 
C CP 
C PDATM 
C PAB1 
C 

C PAB2 
C 

C CONL 
C 

C CONU 
C 

C NPROD 
C MXPRD 

C MXSLAB-  MAXIMUM  NUMBER  OF  SLABS  CURRENTLY  AVAILABLE 
C EPSP  - ERROR  TOLERANCE  FOR  PRESSURES  AT  FLOOR 

C 


- HEIGHT  OF  FLOOR  ABOVE  DATUM  ELEVATION  [M] 

- HEIGHT  OF  LAYER  ABOVE  DATUM  ELEVATION  [M] 

- UPPER  LAYER  TEMPERATURE  [K] 

- LOWER  LAYER  TEMPERATURE  [K] 

- LOWER  LAYER  DENSITY  [KG/M**3] 

- UPPER  LAYER  DENSITY  [KG/M**3] 

- PRESSURE  AT  FLOOR  ABOVE  DATUM  PRESSURE 

[KG/(M*S**2)  - PASCAL] 

- ELEVATION  OF  TOP  OF  VENT  ABOVE  DATUM  ELEVATION  [M] 

- ELEVATION  OF  BOTTOM  OF  VENT  ABOVE  DATUM  ELEVATION  [M] 

- AREA  OF  THE  VENT  [M**2] 

- VALUES  OF  THE  NELEV  ELEVATIONS  ABOVE  DATUM  ELEVATION  [M] 

- PRESSURE  IN  ROOM  1 - PRESSURE  IN  ROOM  2 AT 

ELEVATIONS  YELEV  [KG/(M*S**2)  - PASCAL] 

- NUMBER  OF  ELEVATIONS  DELINEATING  SLABS 

- SPECIFIC  HEAT  [W*S/(KG*K) ] 

- DATUM  ABSOLUTE  PRESSURE  [KG/(M*S**2)  = PASCAL] 

- ABSOLUTE  PRESSURE  IN  ROOM  1 AT  HEIGHT  YELEV 

[KG/(M*S**2)  - PASCAL] 

- ABSOLUTE  PRESSURE  IN  ROOM  2 AT  HEIGHT  YELEV 

[KG/(M*S**2)  - PASCAL] 

- CONCENTRATION  OF  EACH  PRODUCT  IN  LOWER  LAYER 

[UNIT  OF  PRODUCT/ (KG  LAYER)] 

- CONCENTRATION  OF  EACH  PRODUCT  IN  UPPER  LAYER 

[UNIT  OF  PRODUCT/ (KG  LAYER)] 

- NUMBER  OF  PRODUCTS  IN  CURRENT  SCENARIO 

- MAXIMUM  NUMBER  OF  PRODUCTS  CURRENTLY  AVAILABLE 
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C OUTPUT 

C 

C CS1AB 
C 

C PS  LAB 
C 

C PVAB1 
C 

C PVAB2 
C 

C QSLAB 
C DIRS12 
C 

C DPV1M2 
C 

C RSLAB 
C TSLAB 
C YSIAB 
C 

C YVELEV 
C 

C XMSLAB 
C NVELEV 
C NS  LAB 

C 


CONCENTRATION  OF  OTHER  PRODUCTS  IN  EACH  SLAB 
[UNIT  PRODUCT/ (KG  SLAB)] 

AMOUNT  OF  OTHER  PRODUCTS  IN  EACH  SLAB  [UNIT  OF 
PRODUCT/S ] 

ABSOLUTE  PRESSURE  IN  ROOM  1 AT  ELEVATION  YVELEV 
[KG/(M*S**2)  = PASCAL] 

ABSOLUTE  PRESSURE  IN  ROOM  2 AT  ELEVATION  YVELEV 
[KG/(M*S**2)  - PASCAL] 

ENTHALPY  FLOW  RATE  IN  EACH  SLAB  [W] 

A MEASURE  OF  THE  DIRECTION  OF  THE  ROOM  1 TO  ROOM 
2 FLOW  IN  EACH  SLAB 

PRESSURE  IN  ROOM  1 - PRESSURE  IN  ROOM  2 AT 

ELEVATIONS  YVELEV  [KG/(M*S**2)  = PASCAL] 
DENSITY  OF  THE  FLOW  IN  EACH  SLAB  [KG/M**3] 

ABSOLUTE  TEMPERATURE  OF  THE  FLOW  IN  EACH  SLAB  [K] 
ELEVATIONS  ABOVE  THE  DATUM  ELEVATION  OF  THE 

CENTROIDS  OF  MOMENTUM  OF  EACH  SLAB  [M] 
ELEVATIONS  ABOVE  THE  DATUM  ELEVATIONS  OF  VENT 

BOUNDARIES,  LAYERS,  AND  NEUTRAL  PLANES  [M] 
MAGNITUDE  OF  THE  MASS  FLOW  RATE  IN  SLABS  [KG/S] 

NUMBER  OF  UNIQUE  ELEVATIONS  DELINEATING  SLABS 
NUMBER  OF  SLABS  BETWEEN  BOTTOM  AND  TOP  OF  THE  VENT 


C***  ROUTINES  CALLED 
C DELP , HSORT , RMVDUP 

C 
C 

C*END 


DIMENSION  YFLOR(*) , YLAY(*) , TU(*) , TL(*) , DENL(*) , DENU(*) 
DIMENSION  PAB1(*),  PAB2(*),  PFLOR(*) 

DIMENSION  PSLAB(MXSLAB,*) , PV1(*),  PV2(*) 

DIMENSION  CSLAB(MXSLAB,*) , CONL(MXPRD , 2) , CONU(MXPRD , 2) 

DIMENSION  YELEV (*) , DP1M2(*),  DPV1M2(*),  XMSLAB(*) 

DIMENSION  RSLAB(*),  TSLAB (*) , YSLAB(*) , YVELEV (*) 

DIMENSION  QSLAB(*) 

DIMENSION  VP1RT (4) ,VP2RT(4) 

INTEGER  DIRS12(*) 

CHARACTER*]. 32  MSSG 
PARAMETER  ( GAM-1. 40D0) 

DATA  GAMCUT , GAMMAX/ . 471718212, .684731456/ 

C 

NVELEV  - 0 
C 

C***  COPY  VALUES  OF  YELEV(I)  INTO  YV  THAT  LIE  BETWEEN  THE  BOTTOM  AND  TOP 
C OF  THE  VENT 
C 


YMIN  - MAX ( YVBOT , YELEV ( 1 ) ) 
YMAX  - MIN (YVTOP, YELEV (NELEV) ) 
DO  10  I - 1,  NELEV 
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IF (YMIN . LE . YELEV ( I ) . AND . YELEV ( I ) . LE.YMAX)THEN 
NVELEV  - NVELEV  + 1 
YVELEV ( NVELEV ) - YELEV ( I ) 

DPV1M2 (NVELEV)  - DP1M2(I) 

END  IF 
10  CONTINUE 
C 

C***  PUT  YVTOP  INTO  YVELEV  IF  IT  LIES  BETWEEN  YMIN  AND  YMAX 
C AND  CALCULATE  THE  PRESSURE  AT  YVTOP 
C 

IF (YMIN . LE . YVBOT . AND . YVBOT . LE . YMAX) THEN 
NVELEV  - NVELEV  + 1 
YVELEV (NVELEV)  - YVBOT 
CALL  DELP( 

I YVELEV (NVELEV) , YFLOR , YLAY , DENL , DENU , PFLOR , PDATM , 

0 DPV1M2 (NVELEV) , PAB1 (NVELEV) , PAB2 (NVELEV) ) 

END  IF 

C 

C***  PUT  YVBOT  INTO  YVELEV  IF  IT  LIES  BETWEEN  YMIN  AND  YMAX 
C AND  CALCULATE  THE  PRESSURE  AT  YVBOT 
C 

I F (YMIN . LE . YVTOP . AND . YVTOP . LE . YMAX ) THEN 
NVELEV  - NVELEV  + 1 
YVELEV (NVELEV)  - YVTOP 
CALL  DELP( 

1 YVELEV (NVELEV) , YFLOR , YLAY , DENL , DENU , PFLOR , PDATM , 

0 DPV1M2 (NVELEV) , PAB1 (NVELEV) , PAB2 (NVELEV) ) 

ENDIF 

C 

C***  SORT  THE  LIST  AND  REMOVE  DUPLICATES 
C 

CALL  HSORT (YVELEV , DPV1M2 , PABl , PAB2 , NVELEV) 

CALL  RMVDUP (YVELEV , DPV1M2 , PABl , PAB2 , NVELEV) 

C 

C***  THE  NUMBER  OF  SLABS  IS  1 LESS  THAN  THE  NUMBER  OF  HEIGHTS 

C 

NS LAB  - NVELEV  - 1 
DO  30  N - 1,  NSLAB 
C 

C***  DETERMINE  WHETHER  TEMPERATURE  AND  DENSITY  PROPERTIES 
C SHOULD  COME  FROM  ROOM  1 OR  ROOM  2 
C 

PTEST  - DPV1M2(N+1)  + DPV1M2(N) 

IF(PTEST.GT. . 0 ) THEN 
JROOM  - 1 
DIRS12 (N)  - 1 
PRAT  - PAB2 (N) /PABl (N) 

ELSE  IF (PTEST.LT. 0.0) THEN 
DIRS12 (N)  - -1 
JROOM  - 2 

PRAT  - PABl (N) /PAB2 (N) 

ELSE 
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DIRS12 (N)  - 0 
JROOM  - 1 
PRAT  - 0.0D0 
END  IF 
C 

C***  DETERMINE  WHETHER  TEMPERATURE  AND  DENSITY  PROPERTIES 
C SHOULD  COME  FROM  UPPER  OR  LOWER  LAYER 
C 

YMID  = (YVELEV(N)+YVELEV(N+1) )/2 . 0D0 
IF ( YMID . LE . YLAY (JROOM) ) THEN 
TSLAB(N)  - TL( JROOM) 

RSLAB(N)  - DENL ( JROOM )*PRAT 
DO  200  IPROD  - 1,  NPROD 

CSLAB(N, IPROD)  - CONL( IPROD .JROOM) 

200  CONTINUE 

ELSE 

TSLAB(N)  = TU (JROOM) 

RSLAB(N)  = DENU (JROOM) *PRAT 
DO  210  IPROD  = 1,  NPROD 

CSLAB(N, IPROD)  = CONU( IPROD, JROOM) 

210  CONTINUE 

END  IF 

IF(RSLAB(N) . LT . 0 . ) THEN 
WRITE (MSSG , 220)RSLAB(N) 

220  FORMAT (' IN  VENTHP:  RSLAB=' , Ell . 4 , '<0 . ' ) 

CALL  MSGPRT (MSSG , 3 ) 

ENDIF 

C 

C***  FOR  NONZERO -FLOW  SLABS  DETERMINE  XMSLAB(N)  AND  YSLAB(N) 

C 

XMSLAB(N)  = 0.0D0 
QSLAB(N)  - 0.D0 
DO  11  IPROD  - 1,  NPROD 
PSLAB(N, IPROD)  - 0 . 0D0 
11  CONTINUE 

YSLAB(N)  - YMID 
PI  - ABS (DPV1M2 (N) ) 

P2  - ABS(DPV1M2(N+1) ) 

P1RT  - SQRT(Pl) 

P2RT  - SQRT(P2) 

C 

C***  IF  BOTH  CROSS  PRESSURES  ARE  0 THEN  THEN  THERE  IS  NO  FLOW 
C 

IF(P1 . GT . 0 . OR. P2 . GT . 0 . )THEN 
R1  - RSLAB(N) 

Y2  - YVELEV (N+l) 

Y1  - YVELEV (N) 

PABST  - PAB1 (N+l) 

IF(PAB2(N+1) .GT. PABST)PABST  - PAB2(N+1) 

PABSB  - PAB1(N) 

IF(PAB2(N) .GT.PABSB)PABSB  - PAB2(N) 

EPS  - (P1+P2)/(PABST+PABSB) 
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XX  - l.ODO  - EPS 
CC  - . 68DO+. 17DO*EPS 
XO  - l.ODO 

EPSCUT  - EPSP*MAX(XO , PFLOR(l) , PFLOR(2) ) 

EPSCUT  - SQRT( EPSCUT) 

PSUM  - P1RT  + P2RT 
IF(PSUM . EQ . 0 . ODO)THEN 
WP  = O.ODO 
ELSE 

WP  = (PI  + P1RT*P2RT  + P2)/PSUM 
END  IF 
ZZ  - l.ODO 

IF(ABS (WP/EPSCUT) .LE.130)ZZ  = l.DO  - EXP( -WP/EPSCUT) 

C 

C***  GAMCUT  - (2 . /(GAM+1) )**(GAM/(GAM-1 . ) ) 

GAMMAX  - SQRT(  GAM*(2 ./(GAM+1) )**( (GAM+1)/(GAM-1) ) ) 

IF (EPS . LE . 1 . E- 5) THEN 

WW  = l.ODO  - .75DO*EPS/GAM 
ELSE 

IF( EPS. LT. GAMCUT) THEN 
GG  = XX** (l.ODO /GAM) 

FF  - SQRT(  (2 . ODO*GAM/(GAM- 1 . ODO) )*GG*GG* 

$ (1 . ODO-XX/GG)  ) 

ELSE 

FF  - GAMMAX 
END  IF 

WW  - FF/SQRT ( EPS+EPS ) 

ENDIF 

XMSLAB(N)  - CC*WW*SQRT(8*R1)*AVENT*(Y2-Y1) /(YMAX-YMIN)* 

1 ( P2+P1RT*P2RT+P1 ) / ( P2RT  + P1RT) 

XMSLAB(N)  - XMSLAB(N)*ZZ 
XMSLAB(N)  - XMSLAB(N) /3 . ODO 
QSLAB(N)  - CP*XMSLAB(N)*TSLAB(N) 

DO  12  IPROD  - 1,  NPROD 

PSLAB(N, IPROD)  - CSLAB(N, IPROD) *XMSLAB(N) 

12  CONTINUE 

C 

C***  IF  PRESSURES  ARE  CLOSE  TOGETHER  THEN  SET  SLAB  HEIGHT  AT  MID  POINT 
C 

IF(ABS (P1-P2) .GT. . 001) THEN 
C 

C A** 5 - B**5  A**4  + A**3*B  + A**2*B**2  + A*B**3  + B**4 

C - 

C A** 3 - B**3  A** 2 + A*B  + B**2 

C 

C 

C WHEN  A AND  B ARE  "ALMOST"  THE  SAME  THE  RHS  ABOVE  IS  A BETTER  WAY 

C OF  PERFORMING  THE  DESIRED  CALCULATION  SINCE  IT  AVOIDS  CANCELLATION. 

C 

VP1RT(1)  - P1RT 
VP2RT(1)  - P2RT 
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DO  100  I - 2,  4 

VPIRT(I)  - VP1RT ( I - 1 ) *P1RT 
VP2RT(I)  - VP2RT ( I - 1 ) *P2RT 
100  CONTINUE 

RNUM  - VP2RT(4)  + VP1RT(4) 

DO  110  I - 1,  3 

RNUM  - RNUM  + VP1RT(I)*VP2RT(4-I) 

110  CONTINUE 

FACT1  - RNUM/ (PI  + P1RT*P2RT  + P2) 

YSLAB(N)  = ( . 6D0*FACT1*(Y2-Y1)+(P2*Y1-P1*Y2) )/(P2-Pl) 
ELSE 

YSLAB(N)  - YMID 
END  IF 
END  IF 
30  CONTINUE 
RETURN 
END 
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APPENDIX  B - Utility  Algorithms/Subroutines 

The  subroutines  in  this  appendix  are  utility  routines  used  by  some  of  the 
algorithms  presented  in  Appendix  A.  Their  usage  is  explained  in  the  comments 
found  at  the  beginning  of  each  subroutine.  These  routines  do  not  calculate 
physical  quantites  but  manipulate  data  structures  used  by  physical  routines. 
For  example,  HSORT,  sorts  arrays  used  by  the  subroutines  C0MWL1  and  VENTHP. 


. 


HSORT  - Arranges  an  Array  in  Ascending  Order  Using  a Heap  Sorting  Algorithm 


SUBROUTINE  HSORT (RA , RB , RC , RD , N) 

IMPLICIT  DOUBLE  PRECISION  (A-H.O-Z) 

C 

C*BEG 

C***  HSORT  GENERIC  - THIS  ROUTINE  SORTS  THE  ARRAY,  RA,  IN  ASCENDING 
C ORDER  USING  A HEAP  SORT.  CORRESPONDING  CHANGES  ARE  MADE 

C IN  THE  ARRAYS  RB,  RC  AND  RD.  THIS  SUBROUTINE  IS  BASED  ON 

C A SUBROUTINE  NAMED  SORT 2 FOUND  IN: 

C 

C NUMERICAL  RECIPES,  THE  ART  OF  SCIENTIFIC  COMPUTING,  PRESS, 

C WILLIAM  H.,  FLANNERY,  BRAIN  P.,  TEUKOLSKY,  SAUL  AND 

C VETTERLING,  WILLIAM  T. , CAMBRIDGE,  PP.  231-232,  1986 

C 

C*END 

C 

DIMENSION  RA(*) ,RB(*) ,RC(*) ,RD(*) 

L = N/2  + 1 
IR  - N 
10  CONTINUE 

IF(L.GT. 1)THEN 
L = L - 1 
RRA  - RA(L) 

RRB  - RB(L) 

RRC  = RC(L) 

RRD  - RD(L) 

ELSE 

RRA  - RA(IR) 

RRB  - RB(IR) 

RA(IR)  - RA( 1) 

RB(IR)  - RB(1) 

RRC  - RC(IR) 

RRD  - RD(IR) 

RC(IR)  - RC (1) 

RD(IR)  « RD(1) 

IR  - IR  - 1 
IF(IR.EQ.1)THEN 
RA( 1)  - RRA 
RB(1)  - RRB 
RC(1)  - RRC 
RD(1)  - RRD 
RETURN 
END  IF 
END  IF 
I = L 
J - L + L 

20  IF ( J . LE . IR) THEN 

IF(J.LT.IR)THEN 

IF(RA(J).LT.RA(J+1))J  - J + 1 
END  IF 

- 1 


HSORT 


IF(RRA.LT.RA(J) )THEN 
RA( I ) - RA(J) 
RB(I)  = RB(J) 
RC(I)  = RC(J) 
RD(I)  - RD(J) 

I ” J 
J - J + J 
ELSE 

J - IR  + 1 
END  IF 
GO  TO  20 
END  IF 

RA( I ) - RRA 
RB(I)  - RRB 
RC(I)  - RRC 
RD(I)  - RRD 
GO  TO  10 

END 
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INTERP  - Interpolates  A Table  of  Numbers 


SUBROUTINE  INTERP (X , Y , N , T , YINT , YDERIV , I CODE ) 
IMPLICIT  DOUBLE  PRECISION  (A-H.O-Z) 

C*BEG 

C 


C***  INTERP 
C 
C 

C***  SUBROUTINE  ARGUMENTS 
C 


GENERIC  - THIS  ROUTINE  INTERPOLATES  A TABLE  OF 
NUMBERS  FOUND  IN  THE  ARRAYS,  X AND  Y. 


C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 

C*END 


INPUT 


X,  Y 
ICODE 


OUTPUT 


YINT 

YDERIV 


ARRAYS  OF  SIZE  N TO  BE  INTERPOLATED  AT  X-T 
CODE  TO  SELECT  HOW  TO  EXTRAPOLATE  VALUES  IF  T 
IS  LESS  THAN  X(l)  OR  GREATER  THAN  X(N) . 

IF  ICODE  = 1 THEN  YINT  = Y(l)  FOR  T < X(l) 

AND  YINT  - Y(N)  FOR  T > X(N) . 

IF  ICODE  = 2 THEN  YINT  IS  EVALUATED  BY  INTERPOLATION 
IF  X(l)  < T < X(N)  AND  BY  EXTRAPOLATION  IF  T < X(l) 
OR  T > X(N) 


INTERPOLATED  VALUE  OF  THE  Y ARRAY  AT  T 
DERIVATIVE  OF  Y AT  T.  (SLOPE  OF  LINE  SEGMENT 
THAT  BRACKETS  T 


DIMENSION  X(*) , Y(*) 
SAVE 

DATA  I LAST  /l/ 

IF(N. EQ. 1)THEN 
YINT  - Y(l) 

YDERIV  = 0.0D0 
RETURN 
END  IF 

IF(T.LE.X(1))THEN 
IF( ICODE. EQ. 1)THEN 
YINT  - Y(l) 
YDERIV  - 0.0D0 
RETURN 
ELSE 

IMID  - 1 
GO  TO  100 
END  IF 
ENDIF 

IF(T.GE.X(N) )THEN 
IF ( ICODE . EQ . 1 ) THEN 
YINT  - Y(N) 
YDERIV  - 0.0D0 
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RETURN 

ELSE 

IMID  - N - 1 
GO  TO  100 
END  IF 
END  IF 

IF ( ILAST+1 . LE . N) THEN 
IMID  - I LAST 

IF (X(IMID) . LE . T . AND . T . LE .X(IMID+1) )GO  TO  100 
END  IF 

IF ( ILAST+2 . LE . N) THEN 
IMID  - HAST  + 1 

IF(X(IMID) . LE . T . AND . T . LE .X(IMID+1) )GO  TO  100 
END  IF 
IA  - 1 
IZ  - N-l 
1 CONTINUE 

IMID  - (IA+IZ) /2 
IF(T.LT.X( IMID)) THEN 
IZ  - IMID  - 1 
GO  TO  1 
END  IF 

IF ( T . GE . X ( IMID+1 ) ) THEN 
IA  - IMID  + 1 
GO  TO  1 
ENDIF 

100  CONTINUE 

DYDX  - (Y(IMID+1) -Y(IMID) )/(X(IMID+l) -X(IMID) ) 
YINT  - Y(IMID)  + DYDX* (T-X( IMID)) 

YDERIV  - DYDX 
HAST  - IMID 
RETURN 
END 
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RMVDUP  - Remove  Duplicate  Entries  From  a List  of  Data 


SUBROUTINE  RMVDUP (X, Y, Z.W.N) 

IMPLICIT  DOUBLE  PRECISION  (A-H.O-Z) 

RMVDUP  SPECIFIC  - THIS  ROUTINE  REMOVES  DUPLICATE  ENTRIES  FROM 

THE  SORTED  DATA  QUAD-RUPPLES  (X(I) ,Y(I) ,Z(I) ,W(I)  1=1  ...  N) 
THE  PARAMETER  N IS  THEN  ADJUSTED  TO  THE  NEW  NUMBER  OF  ENTRIES. 
TWO  ENTRIES  I AND  1+1  ARE  DUPLICATES  IF  X(I)=X(I+1) 


DIMENSION  X(*),  Y(*),  Z(*)  , W(*) 
II  = 1 

DO  10  I = 2,  N 
DIFF  = X(I)  - X(II) 

IF(DIFF. LT . 0 . )DIFF  = -DIFF 
IF(DIFF.NE.O. )THEN 
II  = II  + 1 
X(II)  = X(I) 

Y(II)  - Y(I) 

Z(II)  = Z(I) 

W(II)  - W(I) 

END  IF 
10  CONTINUE 
N - II 
RETURN 
END 


C*BEG 

C*** 


C 

C*END 
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RMVDUP 


Remove  Duplicate  Entries  From  a List  of  Data 


SUBROUTINE  RMVDUP (X , Y , Z , W , N) 

C 

IMPLICIT  DOUBLE  PRECISION  (A-H.O-Z) 

C 

C*BEG 

C ***  RMVDUP  SPECIFIC  - THIS  ROUTINE  REMOVES  DUPLICATE  ENTRIES  FROM 
C THE  SORTED  DATA  QUAD -RUPP LES  (X(I) , Y(I) , Z(I) ,W(I)  1=1  ...  N) 

C THE  PARAMETER  N IS  THEN  ADJUSTED  TO  THE  NEW  NUMBER  OF  ENTRIES. 

C TWO  ENTRIES  I AND  1+1  ARE  DUPLICATES  IF  X(I)=X(I+1) 

C 

C*END 

C 

C***  LOCAL/DUMMY  VARIABLE  DECLARATIONS 

C 

DIMENSION  X(*),  Y(*)f  Z(*) , W(*) 

II  - 1 

DO  10  I = 2,  N 
DIFF  = X(I)  - X(II) 

IF(DIFF. LT . 0 . )DIFF  = -DIFF 


IF(DIFF . 

NE. 

. 0 . ) THEN 

II  = 

II 

+ 1 

X(II) 

= 

X(I) 

Y(II) 

= 

Y(I) 

Z(II) 

= 

Z(I) 

W(II) 

= 

W(I) 

END  IF 


10  CONTINUE 
N = II 
RETURN 
END 
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